http://wiki.paparazziuav.org/w/api.php?action=feedcontributions&user=Hjkarssies&feedformat=atom
PaparazziUAV - User contributions [en]
2024-03-28T13:13:24Z
User contributions
MediaWiki 1.37.1
http://wiki.paparazziuav.org/w/index.php?title=Working_with_INDI&diff=25217
Working with INDI
2020-10-24T22:24:42Z
<p>Hjkarssies: /* Actuator Dynamics */</p>
<hr />
<div>== Example Airframe Section == <br />
<br />
To use INDI for your drone, you have to specify:<br />
<module name="stabilization" type="indi_simple"/><br />
or<br />
<module name="stabilization" type="rate_indi"/><br />
if you want to fly manual rate control.<br />
<br />
Below you can see an example of how INDI can be configured in the airframe file:<br />
{{Box Code|conf/airframes/myplane.xml|<br />
<source lang="xml"><br />
<section name="STABILIZATION_ATTITUDE_INDI" prefix="STABILIZATION_INDI_"><br />
<!-- control effectiveness --><br />
<define name="G1_P" value="0.0639"/><br />
<define name="G1_Q" value="0.0361"/><br />
<define name="G1_R" value="0.0022"/><br />
<define name="G2_R" value="0.1450"/><br />
<br />
<!-- reference acceleration for attitude control --><br />
<define name="REF_ERR_P" value="600.0"/><br />
<define name="REF_ERR_Q" value="600.0"/><br />
<define name="REF_ERR_R" value="600.0"/><br />
<define name="REF_RATE_P" value="28.0"/><br />
<define name="REF_RATE_Q" value="28.0"/><br />
<define name="REF_RATE_R" value="28.0"/><br />
<br />
<!--Maxium yaw rate, to avoid instability--><br />
<define name="MAX_R" value="120.0" unit="deg/s"/><br />
<br />
<!-- second order filter parameters --><br />
<define name="FILT_CUTOFF" value="8.0"/><br />
<define name="ESTIMATION_FILT_CUTOFF" value="8.0"/><br />
<br />
<!-- first order actuator dynamics --><br />
<define name="ACT_DYN_P" value="0.1"/><br />
<define name="ACT_DYN_Q" value="0.1"/><br />
<define name="ACT_DYN_R" value="0.1"/><br />
<br />
<!-- Adaptive Learning Rate --><br />
<define name="USE_ADAPTIVE" value="FALSE"/><br />
<define name="ADAPTIVE_MU" value="0.0001"/><br />
</section><br />
</source><br />
}}<br />
<br />
Make sure that in your configuration (in paparazzi center) the settings for different controllers, that are not used, are ''removed''. The variables in these settings are not declared, so they will cause errors.<br />
<br />
== Control Effectiveness ==<br />
<br />
The control effectiveness (G1) is what relates inputs to the actuators to angular accelerations. These values are not meant to be tuned. They can be obtained from flight logs, or by using the adaptive algorithm. Note that if you change the inertia of the drone, by adding weight, the control effectiveness is probably different than before. This is also the case when you add or remove bumpers/a flight hull.<br />
<br />
Also note that for quadrotors, a G2 value is introduced, which is there to account for the inertia of the rotors. When the rotors have a certain RPM, they have a certain drag, which is why they provide a yawing moment. More RPM, more drag, more moment. But to change the RPM of the rotors, the motors need to apply a moment to accelerate the rotors. This moment is applied to the drone as well (action = -reaction). To obtain best performance, we have to account for this effect in the controller. That is why there is a control effectiveness value for it.<br />
<br />
== Attitude Control ==<br />
The INDI controller itself controls angular accelerations. This means that if we want to command angular rates or attitude angles, we need to come up with a angular acceleration reference (or virtual control) that will result in the desired rate/attitude. This can be done with a simple PD controller. Since the control effectiveness (which includes the inertia of the vehicle) is accounted for by INDI, the P and D gains do not depend on the size of the drone. Instead, they only depend on the actuator dynamics. <br />
<br />
== Maximum Yaw Rate ==<br />
The yaw control is usually least effective on a quadrotor. Therefore, the yaw command can easily saturate. To avoid yaw overshoot because of this, you can specify a maximum yaw rate.<br />
<br />
== Filtering ==<br />
The angular acceleration is derived from the angular rate. In taking the derivative, noise is amplified. Therefore, the angular acceleration signal can be quite noisy, and filtering is practically always needed. For this purpose a second order filter is used. Such a filter adds a bit of delay to the signal. More filtering will give more delay. This delay will cause disturbances to be measured slightly later and hence be compensated later. If you are not logging the flight data and able to inspect what level of filtering is necessary, it is best not to change this.<br />
<br />
== Actuator Dynamics ==<br />
<br />
Actuators don't react instantly. Motors need time to spin up, and servos can only move so many degrees per second. This means that if the input is suddenly changed, the angular acceleration is not immediately expected to change. If there is actuator feedback available, we know exactly how much the actuator has moved and at what moment how much angular acceleration we should expect. If actuator feedback is not available, an estimate of the actuator position has to be made.<br />
<br />
For electric motors, a first order system can be an accurate model. This can be written as: act(i+1) = act(i) + alpha*(input - act(i)), where alpha is the value that can be specified in the airframe file. Note that the value of alpha depends on the frequency that this model runs! There are a few vehicles for which this value has been determined, such as the ARDrone, Bebop and AscTec Hummingbird. It can be obtained with a step response of the motors, using for instance a logic analyzer and a light sensor to record the spinup profile.<br />
<br />
For more information about estimating actuator dynamics, take a look at the page [[Actuator Dynamics]].<br />
<br />
== Adaptive INDI ==<br />
'''Always use caution with Adaptive INDI'''<br />
<br />
If you enable adaptive INDI, a Least Mean Squares (LMS) algorithm will estimate the control effectiveness values online. It does this by taking the current estimate of the control effectiveness multiplied with the most recent change in input. This is the expected change in angular acceleration. This is then compared to the actual angular acceleration, and the control effectiveness is changed proportional to the error. The parameter ADAPTIVE_MU is the learning rate of the algorithm. A high value will adapt quicker, but can also respond more to disturbances and noise.<br />
<br />
Note that this assumes that all angular accelerations are caused by the inputs. This assumption usually holds indoors, but outside the wind can cause part of the angular acceleration. This can result in faulty adaptation. Also, when taking off with adaptation enabled, faulty learning can occur. This is because while on the ground, inputs will not result in angular accelerations. The algorithm may then estimate a very low control effectiveness. Lastly, if you have no actuator feedback/a bad estimate of the actuator dynamics, the learning may deteriorate.<br />
<br />
One approach is to only use the adaptation for estimation of the control effectiveness after airframe changes. Start flying the drone, and activate the learning while flying through the INDI settings. More importantly, de-activate the learning before landing, or kill the drone before touching the ground. Again, interaction with the ground may result in faulty learning. Save the control effectiveness values.<br />
<br />
== My drone is oscillating ==<br />
If your drone is oscillating or unstable, make sure you check the following (in this order):<br />
<br />
=== Changed inertia ===<br />
If you were flying before, but added something to the drone and now it is oscillating, you probably just need a different control effectiveness. You can try the online adaptation. Make sure that things are not loosely attached to the drone, as this will introduce extra dynamics. Rigid connections are preferred.<br />
<br />
=== The Actuator dynamics ===<br />
If you do not have actuator feedback, it is important to have an accurate actuator model. If you are unsure about your actuator model, it can be worthwhile to try to determine it.<br />
<br />
=== INDI loop vs PD controller ===<br />
<br />
Check if the oscillation comes from the INDI loop (control of angular acceleration) or from the PD controller that commands the reference angular acceleration. It can be that the angular acceleration reference is just too demanding.<br />
<br />
Do this by reducing both the P and D gains. If this does not resolve the issue, you can try holding the drone in your hand (if this can be done safe!) and reducing both P and D gains to zero. The drone will now try to keep zero angular acceleration. If the drone is still oscillating, probably the control effectiveness is wrong (see step 1). Else, try to find P and D values that give a fast response without oscillation.<br />
<br />
== 'Full' INDI ==<br />
<br />
The control described above is the slightly simplified version of INDI (indi_simple), because the controller can only give commands (roll, pitch, yaw, thrust), and does not have full control over each actuator.<br />
The 'control allocation' (which actuators to use for each controlled axis) is done with the static 'motor mixing matrix'.<br />
This is fine in most cases, but it can sometimes lead to problems.<br />
<br />
=== Limitations of indi_simple ===<br />
<br />
First and foremost, the effects of actuator saturation can be unpredictable.<br />
The controller outputs commands, but if some actuators saturate, these commands may not be achievable.<br />
Some commands have more priority than others, pitch and roll are usually more important than yaw for example.<br />
Because the actuators and the controller are decoupled in this setup, dealing with saturation is not possible.<br />
<br />
Second, if the control effectiveness of some of the actuators changes during flight, their new control effectiveness can be estimated with online parameter estimation.<br />
To use this new control effectiveness, there must be an online control allocation, as opposed to a static motor mixing.<br />
<br />
=== Using full INDI ===<br />
<br />
To give INDI access to the individual actuators, the command laws for the servos should use the array <code>actuators_pprz</code>, instead of the motor mixing.<br />
Example:<br />
<br />
<command_laws><br />
<set servo="TOP_LEFT" value="autopilot_get_motors_on() ? actuators_pprz[0] : -MAX_PPRZ"/><br />
<set servo="TOP_RIGHT" value="autopilot_get_motors_on() ? actuators_pprz[1] : -MAX_PPRZ"/><br />
<set servo="BOTTOM_RIGHT" value="autopilot_get_motors_on() ? actuators_pprz[2] : -MAX_PPRZ"/><br />
<set servo="BOTTOM_LEFT" value="autopilot_get_motors_on() ? actuators_pprz[3] : -MAX_PPRZ"/><br />
</command_laws><br />
<br />
Instead of indi_simple, add the stabilization module indi, with the option of RPM feedback (currently only for the bebop 1 and 2).<br />
<br />
<module name="stabilization" type="indi"><br />
<define name="INDI_RPM_FEEDBACK" value="TRUE"/><br />
</module><br />
<br />
In the STABILIZATION_ATTITUDE_INDI section, now the effect of each motor on each axis needs to be specified, like a matrix.<br />
The order is the same as the 'servo' order.<br />
Actuator dynamics can now be specified for each actuator separately.<br />
<br />
{{Box Code|conf/airframes/myplane.xml|<br />
<source lang="xml"><br />
<section name="STABILIZATION_ATTITUDE_INDI" prefix="STABILIZATION_INDI_"><br />
<!-- control effectiveness --><br />
<define name="G1_ROLL" value="{20 , -20, -20 , 20 }"/><br />
<define name="G1_PITCH" value="{14 , 14, -14 , -14 }"/><br />
<define name="G1_YAW" value="{-1, 1, -1, 1}"/><br />
<define name="G1_THRUST" value="{-.4, -.4, -.4, -.4}"/><br />
<!--Counter torque effect of spinning up a rotor--><br />
<define name="G2" value="{-61.2093, 65.3670, -65.7419, 65.4516}"/><br />
<br />
<!-- reference acceleration for attitude control --><br />
<define name="REF_ERR_P" value="600.0"/><br />
<define name="REF_ERR_Q" value="600.0"/><br />
<define name="REF_ERR_R" value="600.0"/><br />
<define name="REF_RATE_P" value="28.0"/><br />
<define name="REF_RATE_Q" value="28.0"/><br />
<define name="REF_RATE_R" value="28.0"/><br />
<br />
<define name="ESTIMATION_FILT_CUTOFF" value="4.0"/><br />
<define name="FILT_CUTOFF" value="5.0"/><br />
<br />
<!-- first order actuator dynamics --><br />
<define name="ACT_DYN" value="{0.1, 0.1, 0.1, 0.1}"/><br />
<br />
<!-- Adaptive Learning Rate --><br />
<define name="USE_ADAPTIVE" value="TRUE"/><br />
<define name="ADAPTIVE_MU" value="0.0001"/><br />
<br />
<!--Priority for each axis (roll, pitch, yaw and thrust)--><br />
<define name="WLS_PRIORITIES" value="{1000, 1000, 1, 100}"/><br />
</section><br />
</source><br />
}}<br />
The adaptive parameter estimation will estimate each element of the G1 and G2 matrices, and can identify differences in actuator performance that are not obvious visual inspection of the airframe.<br />
<br />
=== WLS control allocation ===<br />
<br />
The control allocation with prioritization is incorporated in the form of an active set algorithm, using weighted least squares to determine priorities.<br />
The priorities can be edited with the <code>WLS_PRIORITIES</code> in the section above.<br />
If you do not define this, the default will prioritize pitch and roll, then thrust, and lastly yaw.<br />
This means that if actuators are saturating because of a yaw moment, the aircraft will gradually let go of the yaw control in order to maintain roll and pitch.<br />
<br />
By default, the WLS control allocation is used.<br />
If you want to use the pseudo-inverse to do the control allocation instead (does not take saturation into account!), then define:<br />
<br />
<define name="STABILIZATION_INDI_ALLOCATION_PSEUDO_INVERSE" value="TRUE"/></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=Actuator_Dynamics&diff=25216
Actuator Dynamics
2020-10-24T22:22:33Z
<p>Hjkarssies: </p>
<hr />
<div>[[File:Rpm_sensor.jpg|An RPM sensor|thumb]]<br />
[[File:Saleae-logic.jpg|A Saleae logic analyzer|thumb]]<br />
[[File:Saleae-analog-example.png|A screenshot from the Saleae logic analyser software|thumb]]<br />
[[File:Actuator_dynamics_matlab.svg|An example of estimated actuator dynamics in Matlab|thumb]]<br />
<br />
Some stabilization modules require an estimation of the current actuator states. Examples are the [[Working with INDI|INDI]] and INCA modules. When actuator feedback is not available, this estimation can be mathematical model of the actuators. For brushless motors, a first order transfer function often suffices as actuator model. This guide describes how to accurately determine a first order approximation for a motor with a propeller. Note that the actuator dynamics vary for each configuration of different motors, propellers and ESCs (Electronic Speed Controllers). Used for this example are:<br />
* [https://www.cobramotorsusa.com/motors-2814-16.html Cobra C-2814/16 brushless motor, KV=1050]<br />
* [https://www.banggood.com/1-Piece-HQProp-7040-7x4-CW-or-CCW-Propeller-Slow-Flyer-2-Blade-Green-Nylon-or-Glass-Fiber-Composite-For-RC-Airplane-Multirotors-p-1550754.html HQProp 7040 7x4 Propeller Slow Flyer 2 Blade Nylon/Glass Fiber Composite]<br />
* [http://kiss.flyduino.net/kiss-esc-24a-re/ KISS ESC 24A re – 32bit Electronic Speed Controller]<br><br><br />
<br />
==Actuator model==<br />
With this method, actuators are approximated by a first order model in de following form:<br><br><br />
<math>H_{act}=\frac{K}{\tau s+1}</math><br><br><br />
where <math>H_{act}</math> is the actuator transfer function in the Laplace domain, <math>K</math> is the system's gain, and <math>\tau</math> is the system's time constant.<br><br><br />
<br />
When using such a model on discrete time systems, an estimation of an actuator's current position (<math>\delta_{est}</math>) can be calculated as a function of the current actuator command (<math>\delta</math>) and the previous estimation of its position (<math>\delta_{prev}</math>):<br><br><br />
<math>\delta_{est}=\delta_{prev}+\alpha(\delta-\delta_{prev})</math><br><br />
where <math>\alpha=1-e^{-\tau\Delta t}</math><br><br />
where <math>\Delta t</math> is the duration of one time step<br><br><br />
<br />
To make a good approximation of actuator dynamics, one needs to determine appropriate values for <math>\tau</math> and <math>\alpha</math><br><br><br />
<br />
==Measuring Actuator Response==<br />
The first step in approximating actuator dynamics is to measure actuator response. This can be done with by connecting an RPM sensor to the motor's wires, and recording both the ESC's input PWM and output RPM signals. An schematic representation of the setup used for this example is shown in the figure below:<br><br><br />
<br />
[[File:Actuator_Dynamics_Measurements_Schematic.svg|A schematic overview of a test setup to measure actuator dynamics|400px]]<br><br><br />
<br />
Used for the setup of this example:<br />
* [https://hobbyking.com/en_us/orange-rpm-sensor.html Orange RPM Sensor] for RPM measurements<br />
* [https://www.saleae.com/ Saleae Logic 8 Digital Analyzer] for recording measurements<br />
* [https://www.mathworks.com/products/matlab.html MATLAB] for analyzing measurements<br><br><br />
<br />
In order to get a representative model for the whole range of the motor, control the motor with varying speeds and varying step sizes in the PWM signal. You can use Paparazzi's module ''servo_tester'', which runs a predefined program. Try to experiment a bit with the module's settings to create a test program suitable for your purpose, and make sure to firmly attach your motor before running the program.<br><br><br />
<br />
You can record the PWM and RPM signals with your digital analyzer. Make sure to select appropriate bitrates. Saleae provides a function in their software with which you can export your measurements to a MATLAB file, which you can read with the MATLAB script below.<br><br><br />
<br />
==Data Processing==<br />
With the following MATLAB script, you can make the estimation of you first order actuator model's time constant. In the script, indicate the file name of your measurements, the channel numbers for the different types of data, and the timeframe in which you recorded your relevant data. You can also specify different values of tau to try. The script will try these values, and interpolate results to find the tau with which the simulated response best fits the measured response. Don't forget to set the correct frequency at which your module runs if you want to get a correct value for alpha as well. The minimal period settings are used to filter out any signal changes that are too quick, and are therefore probably misreadings.<br><br><br />
<br />
<source lang="Matlab"><br />
%{<br />
This script is written for estimation of first order actuator dynamics for<br />
motors with rotors. It uses measurements of both a PWM input signal and the<br />
actuator's RPM response, and tries to fit a first order model to these data<br />
of de form:<br />
<br />
K <br />
-----------<br />
tau * s + 1<br />
<br />
This form can be used to estimate a actuator's state as a function of its<br />
previous state and its current input as follows:<br />
<br />
u_est = u_prev + alpha * (u - u_prev)<br />
where alpha = 1 - e^(tau * dt)<br />
<br />
© H.J. Karssies - Delft University of Technology<br />
%}<br />
<br />
% Close figures and clear workspace<br />
clc<br />
clear all<br />
close all<br />
<br />
% Settings<br />
data_file = 'motor_data_04.mat'; % Data file<br />
pwm.channel_no = 0; % Channel number of PWM input data<br />
pwm.period_min = 0.0005; % Minimum PWM period used for filtering<br />
rpm.channel_no = 1; % Channel number of RPM measurement data<br />
rpm.period_min = 0.00005; % Minimum RPM period used for filtering<br />
t_min = 5; % Starting time of relevant data (s)<br />
t_max = 10.4; % End time of relevant data (s)<br />
taus = 5:5:100; % List of values for tau to try (s)<br />
main_freq = 512; % Paparazzi main frequency<br />
<br />
% Load data<br />
raw_data = load(data_file);<br />
pwm = load_channel(pwm, raw_data);<br />
rpm = load_channel(rpm, raw_data);<br />
pwm.freq = length(pwm.time)/(pwm.time(end)-pwm.time(1));<br />
<br />
% Calculate motor frequencies<br />
pwm.motor_freq = pwm.period*1000000;<br />
rpm.motor_freq = 1./rpm.period;<br />
<br />
% Check performance for different values of tau<br />
for n = 1:length(taus)<br />
errors(n) = calc_performance(taus(n), pwm, rpm, t_min, t_max, false);<br />
end<br />
<br />
% Interpolate results and find minimum<br />
taus_int = taus(1):0.01:taus(end);<br />
errors_int = spline(taus, errors, taus_int);<br />
error_min = min(errors_int);<br />
tau_min = taus_int(find(errors_int == error_min));<br />
alpha = 1-exp(-tau_min/main_freq);<br />
<br />
% Plot result<br />
figure;<br />
set(gcf, 'Position', [100, 100, 1080, 540]);<br />
subplot(1,2,1);<br />
plot(taus_int, errors_int);<br />
hold on<br />
scatter(tau_min, error_min, 'black');<br />
grid on<br />
title("\tau estimation");<br />
xlabel("\tau (s)");<br />
ylabel("Error (rpm²)");<br />
text(tau_min, error_min+0.1*diff(ylim), ...<br />
[sprintf("\\tau = %.2f s", tau_min); ...<br />
sprintf("\\alpha = %.5f (%d Hz)", [alpha, main_freq])]);<br />
<br />
% Perform calculations with final tau and plot figures<br />
calc_performance(tau_min, pwm, rpm, t_min, t_max, true);<br />
<br />
% Load channel data and process it<br />
function channel = load_channel(channel, data)<br />
<br />
% Process channel<br />
channel_data = eval('data.digital_channel_' + string(channel.channel_no));<br />
channel.bit_1 = data.digital_channel_initial_bitstates(channel.channel_no+1);<br />
channel.period = channel_data(2-channel.bit_1:2:end)./data.digital_sample_rate_hz;<br />
channel.time = cumsum(channel_data./data.digital_sample_rate_hz);<br />
channel.time = channel.time(2-channel.bit_1:2:end);<br />
<br />
% Filter out invalid data<br />
channel.time = channel.time(channel.period>channel.period_min);<br />
channel.period = channel.period(channel.period>channel.period_min);<br />
<br />
end<br />
<br />
% Calculate performance with certain tau<br />
function error = calc_performance(tau, pwm, rpm, t_min, t_max, plot_fig)<br />
<br />
% Simulate RPM response<br />
motor = tf(tau,[1 tau]);<br />
<br />
% Run simulation at rounded PWM sample times<br />
bias = mean(pwm.motor_freq(1:10));<br />
dt = 1/round(pwm.freq);<br />
t = [0:length(pwm.time)-1]*dt;<br />
motor = tf(tau,[1 tau]);<br />
sim_pwm_freq = lsim(motor, pwm.motor_freq-bias, t) + bias;<br />
<br />
% Get the RPM requency values at the rounden PWM sample times<br />
int_rpm_freq = interp1(rpm.time,rpm.motor_freq,t,[],1)';<br />
int_rpm_freq(isinf(int_rpm_freq)) = 1;<br />
<br />
% Fit a simple model y = ax + b<br />
data_range = and(t>=t_min, t<=t_max);<br />
inputs = [ones(size(sim_pwm_freq(data_range))) ...<br />
sim_pwm_freq(data_range)];<br />
gain = inputs \ int_rpm_freq(data_range);<br />
sim_rpm_freq = inputs*gain;<br />
<br />
% Calculate squared error<br />
error = sum((int_rpm_freq(data_range)-sim_rpm_freq).^2);<br />
<br />
% Plot figure<br />
if plot_fig<br />
subplot(1,2,2);<br />
plot(pwm.time, pwm.motor_freq)<br />
hold on<br />
plot(rpm.time, rpm.motor_freq)<br />
plot(pwm.time(data_range), sim_rpm_freq)<br />
grid on<br />
xlim([0, max(pwm.time)]);<br />
ylim([0, 1.2*max(rpm.motor_freq)]);<br />
plot([t_min t_min], ylim, 'black');<br />
plot([t_max t_max], ylim, 'black');<br />
title(sprintf("PWM signal and response for \\tau = %.2f s", tau));<br />
xlabel("Time (s)");<br />
ylabel("Duty time (ms) / motor speed (rpm)");<br />
legend(["PWM signal", "RPM measurement", ...<br />
"Response approximation"])<br />
end<br />
<br />
end<br />
</source></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=Actuator_Dynamics&diff=25215
Actuator Dynamics
2020-10-24T21:41:47Z
<p>Hjkarssies: /* Measuring Actuator Response */</p>
<hr />
<div>[[File:Rpm_sensor.jpg|An RPM sensor|thumb]]<br />
[[File:Saleae-logic.jpg|A Saleae logic analyzer|thumb]]<br />
[[File:Saleae-analog-example.png|A screenshot from the Saleae logic analyser software|thumb]]<br />
[[File:Actuator_dynamics_matlab.svg|An example of estimated actuator dynamics in Matlab|thumb]]<br />
<br />
Some stabilization modules require an estimation of the current actuator states. Examples are the [[Working with INDI|INDI]] and INCA modules. When actuator feedback is not available, this estimation can be mathematical model of the actuators. For brushless motors, a first order transfer function often suffices as actuator model. This guide describes how to accurately determine a first order approximation for a motor with a propeller. Note that the actuator dynamics vary for each configuration of different motors, propellers and ESCs (Electronic Speed Controllers). Used for this example are:<br />
* [https://www.cobramotorsusa.com/motors-2814-16.html Cobra C-2814/16 brushless motor, KV=1050]<br />
* [https://www.banggood.com/1-Piece-HQProp-7040-7x4-CW-or-CCW-Propeller-Slow-Flyer-2-Blade-Green-Nylon-or-Glass-Fiber-Composite-For-RC-Airplane-Multirotors-p-1550754.html HQProp 7040 7x4 Propeller Slow Flyer 2 Blade Nylon/Glass Fiber Composite]<br />
* [http://kiss.flyduino.net/kiss-esc-24a-re/ KISS ESC 24A re – 32bit Electronic Speed Controller]<br><br><br />
<br />
==Actuator model==<br />
With this method, actuators are approximated by a first order model in de following form:<br><br><br />
<math>H_{act}=\frac{K}{\tau s+1}</math><br><br><br />
where <math>H_{act}</math> is the actuator transfer function in the Laplace domain, <math>K</math> is the system's gain, and <math>\tau</math> is the system's time constant.<br><br><br />
<br />
When using such a model on discrete time systems, an estimation of an actuator's current position (<math>\delta_{est}</math>) can be calculated as a function of the current actuator command (<math>\delta</math>) and the previous estimation of its position (<math>\delta_{prev}</math>):<br><br><br />
<math>\delta_{est}=\delta_{prev}+\alpha(\delta-\delta_{prev})</math><br><br />
where <math>\alpha=1-e^{-\tau\Delta t}</math><br><br />
where <math>\Delta t</math> is the duration of one time step<br><br><br />
<br />
To make a good approximation of actuator dynamics, one needs to determine appropriate values for <math>\tau</math> and <math>\alpha</math><br><br><br />
<br />
==Measuring Actuator Response==<br />
The first step in approximating actuator dynamics is to measure actuator response. This can be done with by connecting an RPM sensor to the motor's wires, and recording both the ESC's input PWM and output RPM signals. An schematic representation of the setup used for this example is shown in the figure below:<br><br><br />
<br />
[[File:Actuator_Dynamics_Measurements_Schematic.svg|A schematic overview of a test setup to measure actuator dynamics|400px]]<br><br><br />
<br />
Used for the setup of this example:<br />
* [https://hobbyking.com/en_us/orange-rpm-sensor.html Orange RPM Sensor] for RPM measurements<br />
* [https://www.saleae.com/ Saleae Logic 8] for recording measurements<br />
* [https://www.mathworks.com/products/matlab.html MATLAB] for analyzing measurements<br><br><br />
<br />
In order to get a representative model for the whole range of the motor, control the motor with varying speeds and varying step sizes in the PWM signal. You can use Paparazzi's module ''servo_tester'', which runs a predefined program. Try to experiment a bit with the module's settings to create a test program suitable for your purpose, and make sure to firmly attach your motor before running the program.<br />
<br />
==WORK IN PROGRESS==<br />
<br />
<source lang="Matlab"><br />
%{<br />
This script is written for estimation of first order actuator dynamics for<br />
motors with rotors. It uses measurements of both a PWM input signal and the<br />
actuator's RPM response, and tries to fit a first order model to these data<br />
of de form:<br />
<br />
K <br />
-----------<br />
tau * s + 1<br />
<br />
This form can be used to estimate a actuator's state as a function of its<br />
previous state and its current input as follows:<br />
<br />
u_est = u_prev + alpha * (u - u_prev)<br />
where alpha = 1 - e^(tau * dt)<br />
<br />
© H.J. Karssies - Delft University of Technology<br />
%}<br />
<br />
% Close figures and clear workspace<br />
clc<br />
clear all<br />
close all<br />
<br />
% Settings<br />
data_file = 'motor_data_04.mat'; % Data file<br />
pwm.channel_no = 0; % Channel number of PWM input data<br />
pwm.period_min = 0.0005; % Minimum PWM period used for filtering<br />
rpm.channel_no = 1; % Channel number of RPM measurement data<br />
rpm.period_min = 0.00005; % Minimum RPM period used for filtering<br />
t_min = 5; % Starting time of relevant data (s)<br />
t_max = 10.4; % End time of relevant data (s)<br />
taus = 5:5:100; % List of values for tau to try (s)<br />
main_freq = 512; % Paparazzi main frequency<br />
<br />
% Load data<br />
raw_data = load(data_file);<br />
pwm = load_channel(pwm, raw_data);<br />
rpm = load_channel(rpm, raw_data);<br />
pwm.freq = length(pwm.time)/(pwm.time(end)-pwm.time(1));<br />
<br />
% Calculate motor frequencies<br />
pwm.motor_freq = pwm.period*1000000;<br />
rpm.motor_freq = 1./rpm.period;<br />
<br />
% Check performance for different values of tau<br />
for n = 1:length(taus)<br />
errors(n) = calc_performance(taus(n), pwm, rpm, t_min, t_max, false);<br />
end<br />
<br />
% Interpolate results and find minimum<br />
taus_int = taus(1):0.01:taus(end);<br />
errors_int = spline(taus, errors, taus_int);<br />
error_min = min(errors_int);<br />
tau_min = taus_int(find(errors_int == error_min));<br />
alpha = 1-exp(-tau_min/main_freq);<br />
<br />
% Plot result<br />
figure;<br />
set(gcf, 'Position', [100, 100, 1080, 540]);<br />
subplot(1,2,1);<br />
plot(taus_int, errors_int);<br />
hold on<br />
scatter(tau_min, error_min, 'black');<br />
grid on<br />
title("\tau estimation");<br />
xlabel("\tau (s)");<br />
ylabel("Error (rpm²)");<br />
text(tau_min, error_min+0.1*diff(ylim), ...<br />
[sprintf("\\tau = %.2f s", tau_min); ...<br />
sprintf("\\alpha = %.5f (%d Hz)", [alpha, main_freq])]);<br />
<br />
% Perform calculations with final tau and plot figures<br />
calc_performance(tau_min, pwm, rpm, t_min, t_max, true);<br />
<br />
% Load channel data and process it<br />
function channel = load_channel(channel, data)<br />
<br />
% Process channel<br />
channel_data = eval('data.digital_channel_' + string(channel.channel_no));<br />
channel.bit_1 = data.digital_channel_initial_bitstates(channel.channel_no+1);<br />
channel.period = channel_data(2-channel.bit_1:2:end)./data.digital_sample_rate_hz;<br />
channel.time = cumsum(channel_data./data.digital_sample_rate_hz);<br />
channel.time = channel.time(2-channel.bit_1:2:end);<br />
<br />
% Filter out invalid data<br />
channel.time = channel.time(channel.period>channel.period_min);<br />
channel.period = channel.period(channel.period>channel.period_min);<br />
<br />
end<br />
<br />
% Calculate performance with certain tau<br />
function error = calc_performance(tau, pwm, rpm, t_min, t_max, plot_fig)<br />
<br />
% Simulate RPM response<br />
motor = tf(tau,[1 tau]);<br />
<br />
% Run simulation at rounded PWM sample times<br />
bias = mean(pwm.motor_freq(1:10));<br />
dt = 1/round(pwm.freq);<br />
t = [0:length(pwm.time)-1]*dt;<br />
motor = tf(tau,[1 tau]);<br />
sim_pwm_freq = lsim(motor, pwm.motor_freq-bias, t) + bias;<br />
<br />
% Get the RPM requency values at the rounden PWM sample times<br />
int_rpm_freq = interp1(rpm.time,rpm.motor_freq,t,[],1)';<br />
int_rpm_freq(isinf(int_rpm_freq)) = 1;<br />
<br />
% Fit a simple model y = ax + b<br />
data_range = and(t>=t_min, t<=t_max);<br />
inputs = [ones(size(sim_pwm_freq(data_range))) ...<br />
sim_pwm_freq(data_range)];<br />
gain = inputs \ int_rpm_freq(data_range);<br />
sim_rpm_freq = inputs*gain;<br />
<br />
% Calculate squared error<br />
error = sum((int_rpm_freq(data_range)-sim_rpm_freq).^2);<br />
<br />
% Plot figure<br />
if plot_fig<br />
subplot(1,2,2);<br />
plot(pwm.time, pwm.motor_freq)<br />
hold on<br />
plot(rpm.time, rpm.motor_freq)<br />
plot(pwm.time(data_range), sim_rpm_freq)<br />
grid on<br />
xlim([0, max(pwm.time)]);<br />
ylim([0, 1.2*max(rpm.motor_freq)]);<br />
plot([t_min t_min], ylim, 'black');<br />
plot([t_max t_max], ylim, 'black');<br />
title(sprintf("PWM signal and response for \\tau = %.2f s", tau));<br />
xlabel("Time (s)");<br />
ylabel("Duty time (ms) / motor speed (rpm)");<br />
legend(["PWM signal", "RPM measurement", ...<br />
"Response approximation"])<br />
end<br />
<br />
end<br />
</source></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=Iris&diff=25214
Iris
2020-10-24T21:34:53Z
<p>Hjkarssies: /* Calibrate the ESCs */</p>
<hr />
<div>[[File:Iris.jpg|thumb|The 3DR Iris+ with 3DR radio and 3DR RC receiver]]<br />
<br />
<br />
= Introduction =<br />
<br />
The Iris+ drone was manufactured by 3DR and is based around the Pixhawk autopilot. Paparazzi [[Pixhawk|supports]] the Pixhawk hardware, and this page serves as simple step-by-step guide to get a Pixhawk based drone flying on with the Paparazzi autopilot. In this example, only standard Pixhawk hardware is used, as it comes delivered with the stock Iris.<br />
<br />
=Videos=<br />
Paparazzi UAV Iris+ Pixhawk tutorial videos - by the [http://mavlab.lr.tudelft.nl/ TU Delft MAVLab] <br><br />
{{#ev:youtubeplaylist|PL_KSX9GOn2P-CVsy-lRSHAnUwVF9wI_pZ}} <br />
<br style="clear:both"><br />
<br />
= Getting the Iris flight ready=<br />
<br />
You'll need:<br />
# The stock Iris drone, (equipped with the Pixhawk), <br />
# Its stock 3DR radio's and 3DR RC. <br />
# The Pixhawk should still contain the PX4 / APM bootloader. <br />
# A computer with Paparazzi [[Installation|installed]].<br />
# A micro USB cable.<br />
# Optional: the 3DR Tarrot Gimbal and a GoPro 3.<br />
<br />
Steps you will need to perform: <br />
# Plug in the USB cable to the Iris. <b>Remove the props. Leave the battery disconnected.</b><br />
# Flash both the AP and FBW. [[Pixhawk#Flashing|On how to flash, follow this link here]]<br />
# Power cycle the Pixhawk (e.g. replug USB).<br />
# Make sure to set the target to AP again, and build. Uploading is not necessary, this step is required to prevent md5 checksum errors.<br />
# Plug in the 3DR radio to the computer. Select the <i>Flight USB0 @57600</i> session in Paparazzi Center, and click execute. Paparazzi GCS will start.<br />
# [[ImuCalibration|Calibrating the Magnetometer]] and [[#Calibrate the ESCs | calibrate the ESCs]] of the Iris+.<br />
<br />
==Calibrate the ESCs==<br />
The original PX4 firmware calibrates the ESCs automatically, paparazzi does not support this. Therefor it is needed to calibrate the ESCs manually. For this two options are available:<br />
* Calibrate by using a servo tester (see page [[Actuator Dynamics]])<br />
* Calibrate using <i> mode RC direct</i><br />
The first option requires the user to have a servo tester, open up the Iris+ and disconnect wires. Therefor, this wiki describes the second option.<br />
<br />
In order to use the <i>RC direct mode</i>, a change is required in the airframe xml configuration. Comment or temporarily remove the following line:<br />
{{Box Code|conf/airframes/TUDELFT/tudelft_iris_indi.xml|<source lang="xml"><br />
<define name="FBW_MODE_AUTO_ONLY" value="true"/><br />
</source>}}<br />
Optionally: change the following lines in the section AUTOPILOT:<br />
{{Box Code|conf/airframes/TUDELFT/tudelft_iris_indi.xml|<source lang="xml"><br />
<define name="MODE_STARTUP" value="AP_MODE_ATTITUDE_DIRECT" /><br />
<define name="MODE_MANUAL" value="AP_MODE_ATTITUDE_DIRECT" /><br />
</source>}}<br />
to:<br />
{{Box Code|conf/airframes/TUDELFT/tudelft_iris_indi.xml|<source lang="xml"><br />
<define name="MODE_STARTUP" value="AP_MODE_RC_DIRECT" /><br />
<define name="MODE_MANUAL" value="AP_MODE_RC_DIRECT" /><br />
</source>}}<br />
This last step is not required, but properly communicates the new mode to the GCS software.<br />
<br />
<b> Important: for safety it is <i>strongly</i> recommended to remove the propellers </b> <br/><br />
<br />
[[File:3DRRC.jpg|thumb|Close up of the 3DR RC receiver. 1) The manual Kill switch. 2) The Gimbal control. 3) The mode switch.]]<br />
Build and upload the FBW firmware, and if the optional steps were taken, also build and upload the AP. <br />
<br />
The system is now ready to calibrate the ESCs. Perform the following steps:<br />
* Power off the Iris+<br />
* On the RC<br />
** Set the set radio mode (switch 3 in figure 2) to <i>std</i><br />
** Disable manual kill by setting the kill switch (switch 1 in figure 2) to <i>on</i><br />
** Put the throttle to max.<br />
** Turn on the RC<br />
* Plug in the usb cable to the Iris+. The AP can now boot. <br />
* Wait approximately 10 seconds until the leds blinks slow (once per second).<br />
* Plug in the battery to the Iris+ (with props removed!). You should hear a beep coming from the ESCs.<br />
* On the RC, put the throttle to zero. You should hear confirmation beeps from the ESCs.<br />
* Validate the range of the throttle stick by slowly throttling up to max. The motors should spin accordingly.<br />
* Revert the changes made to the airframe xml, and upload the firmware(s).</div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=File:Actuator_Dynamics_Measurements_Schematic.svg&diff=25213
File:Actuator Dynamics Measurements Schematic.svg
2020-10-24T21:31:40Z
<p>Hjkarssies: Hjkarssies uploaded a new version of File:Actuator Dynamics Measurements Schematic.svg</p>
<hr />
<div></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=File:Actuator_Dynamics_Measurements_Schematic.svg&diff=25212
File:Actuator Dynamics Measurements Schematic.svg
2020-10-24T21:30:34Z
<p>Hjkarssies: Hjkarssies uploaded a new version of File:Actuator Dynamics Measurements Schematic.svg</p>
<hr />
<div></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=Actuator_Dynamics&diff=25211
Actuator Dynamics
2020-10-18T20:47:18Z
<p>Hjkarssies: </p>
<hr />
<div>[[File:Rpm_sensor.jpg|An RPM sensor|thumb]]<br />
[[File:Saleae-logic.jpg|A Saleae logic analyzer|thumb]]<br />
[[File:Saleae-analog-example.png|A screenshot from the Saleae logic analyser software|thumb]]<br />
[[File:Actuator_dynamics_matlab.svg|An example of estimated actuator dynamics in Matlab|thumb]]<br />
<br />
Some stabilization modules require an estimation of the current actuator states. Examples are the [[Working with INDI|INDI]] and INCA modules. When actuator feedback is not available, this estimation can be mathematical model of the actuators. For brushless motors, a first order transfer function often suffices as actuator model. This guide describes how to accurately determine a first order approximation for a motor with a propeller. Note that the actuator dynamics vary for each configuration of different motors, propellers and ESCs (Electronic Speed Controllers). Used for this example are:<br />
* [https://www.cobramotorsusa.com/motors-2814-16.html Cobra C-2814/16 brushless motor, KV=1050]<br />
* [https://www.banggood.com/1-Piece-HQProp-7040-7x4-CW-or-CCW-Propeller-Slow-Flyer-2-Blade-Green-Nylon-or-Glass-Fiber-Composite-For-RC-Airplane-Multirotors-p-1550754.html HQProp 7040 7x4 Propeller Slow Flyer 2 Blade Nylon/Glass Fiber Composite]<br />
* [http://kiss.flyduino.net/kiss-esc-24a-re/ KISS ESC 24A re – 32bit Electronic Speed Controller]<br><br><br />
<br />
==Actuator model==<br />
With this method, actuators are approximated by a first order model in de following form:<br><br><br />
<math>H_{act}=\frac{K}{\tau s+1}</math><br><br><br />
where <math>H_{act}</math> is the actuator transfer function in the Laplace domain, <math>K</math> is the system's gain, and <math>\tau</math> is the system's time constant.<br><br><br />
<br />
When using such a model on discrete time systems, an estimation of an actuator's current position (<math>\delta_{est}</math>) can be calculated as a function of the current actuator command (<math>\delta</math>) and the previous estimation of its position (<math>\delta_{prev}</math>):<br><br><br />
<math>\delta_{est}=\delta_{prev}+\alpha(\delta-\delta_{prev})</math><br><br />
where <math>\alpha=1-e^{-\tau\Delta t}</math><br><br />
where <math>\Delta t</math> is the duration of one time step<br><br><br />
<br />
To make a good approximation of actuator dynamics, one needs to determine appropriate values for <math>\tau</math> and <math>\alpha</math><br><br><br />
<br />
==Measuring Actuator Response==<br />
The first step in approximating actuator dynamics is to measure actuator response. This can be done with by connecting an RPM sensor to the motor's wires, and recording both the ESC's input PWM and output RPM signals. An schematic representation of the setup used for this example is shown in the figure below:<br><br><br />
<br />
[[File:Actuator_Dynamics_Measurements_Schematic.svg|A schematic overview of a test setup to measure actuator dynamics|400px]]<br><br><br />
<br />
Used for the setup of this example:<br />
* [https://hobbyking.com/en_us/orange-rpm-sensor.html Orange RPM Sensor] for RPM measurements<br />
* [https://www.saleae.com/ Saleae Logic 8] for recording measurements<br />
* [https://www.mathworks.com/products/matlab.html MATLAB] for analyzing measurements<br><br><br />
<br />
==WORK IN PROGRESS==<br />
<br />
<source lang="Matlab"><br />
%{<br />
This script is written for estimation of first order actuator dynamics for<br />
motors with rotors. It uses measurements of both a PWM input signal and the<br />
actuator's RPM response, and tries to fit a first order model to these data<br />
of de form:<br />
<br />
K <br />
-----------<br />
tau * s + 1<br />
<br />
This form can be used to estimate a actuator's state as a function of its<br />
previous state and its current input as follows:<br />
<br />
u_est = u_prev + alpha * (u - u_prev)<br />
where alpha = 1 - e^(tau * dt)<br />
<br />
© H.J. Karssies - Delft University of Technology<br />
%}<br />
<br />
% Close figures and clear workspace<br />
clc<br />
clear all<br />
close all<br />
<br />
% Settings<br />
data_file = 'motor_data_04.mat'; % Data file<br />
pwm.channel_no = 0; % Channel number of PWM input data<br />
pwm.period_min = 0.0005; % Minimum PWM period used for filtering<br />
rpm.channel_no = 1; % Channel number of RPM measurement data<br />
rpm.period_min = 0.00005; % Minimum RPM period used for filtering<br />
t_min = 5; % Starting time of relevant data (s)<br />
t_max = 10.4; % End time of relevant data (s)<br />
taus = 5:5:100; % List of values for tau to try (s)<br />
main_freq = 512; % Paparazzi main frequency<br />
<br />
% Load data<br />
raw_data = load(data_file);<br />
pwm = load_channel(pwm, raw_data);<br />
rpm = load_channel(rpm, raw_data);<br />
pwm.freq = length(pwm.time)/(pwm.time(end)-pwm.time(1));<br />
<br />
% Calculate motor frequencies<br />
pwm.motor_freq = pwm.period*1000000;<br />
rpm.motor_freq = 1./rpm.period;<br />
<br />
% Check performance for different values of tau<br />
for n = 1:length(taus)<br />
errors(n) = calc_performance(taus(n), pwm, rpm, t_min, t_max, false);<br />
end<br />
<br />
% Interpolate results and find minimum<br />
taus_int = taus(1):0.01:taus(end);<br />
errors_int = spline(taus, errors, taus_int);<br />
error_min = min(errors_int);<br />
tau_min = taus_int(find(errors_int == error_min));<br />
alpha = 1-exp(-tau_min/main_freq);<br />
<br />
% Plot result<br />
figure;<br />
set(gcf, 'Position', [100, 100, 1080, 540]);<br />
subplot(1,2,1);<br />
plot(taus_int, errors_int);<br />
hold on<br />
scatter(tau_min, error_min, 'black');<br />
grid on<br />
title("\tau estimation");<br />
xlabel("\tau (s)");<br />
ylabel("Error (rpm²)");<br />
text(tau_min, error_min+0.1*diff(ylim), ...<br />
[sprintf("\\tau = %.2f s", tau_min); ...<br />
sprintf("\\alpha = %.5f (%d Hz)", [alpha, main_freq])]);<br />
<br />
% Perform calculations with final tau and plot figures<br />
calc_performance(tau_min, pwm, rpm, t_min, t_max, true);<br />
<br />
% Load channel data and process it<br />
function channel = load_channel(channel, data)<br />
<br />
% Process channel<br />
channel_data = eval('data.digital_channel_' + string(channel.channel_no));<br />
channel.bit_1 = data.digital_channel_initial_bitstates(channel.channel_no+1);<br />
channel.period = channel_data(2-channel.bit_1:2:end)./data.digital_sample_rate_hz;<br />
channel.time = cumsum(channel_data./data.digital_sample_rate_hz);<br />
channel.time = channel.time(2-channel.bit_1:2:end);<br />
<br />
% Filter out invalid data<br />
channel.time = channel.time(channel.period>channel.period_min);<br />
channel.period = channel.period(channel.period>channel.period_min);<br />
<br />
end<br />
<br />
% Calculate performance with certain tau<br />
function error = calc_performance(tau, pwm, rpm, t_min, t_max, plot_fig)<br />
<br />
% Simulate RPM response<br />
motor = tf(tau,[1 tau]);<br />
<br />
% Run simulation at rounded PWM sample times<br />
bias = mean(pwm.motor_freq(1:10));<br />
dt = 1/round(pwm.freq);<br />
t = [0:length(pwm.time)-1]*dt;<br />
motor = tf(tau,[1 tau]);<br />
sim_pwm_freq = lsim(motor, pwm.motor_freq-bias, t) + bias;<br />
<br />
% Get the RPM requency values at the rounden PWM sample times<br />
int_rpm_freq = interp1(rpm.time,rpm.motor_freq,t,[],1)';<br />
int_rpm_freq(isinf(int_rpm_freq)) = 1;<br />
<br />
% Fit a simple model y = ax + b<br />
data_range = and(t>=t_min, t<=t_max);<br />
inputs = [ones(size(sim_pwm_freq(data_range))) ...<br />
sim_pwm_freq(data_range)];<br />
gain = inputs \ int_rpm_freq(data_range);<br />
sim_rpm_freq = inputs*gain;<br />
<br />
% Calculate squared error<br />
error = sum((int_rpm_freq(data_range)-sim_rpm_freq).^2);<br />
<br />
% Plot figure<br />
if plot_fig<br />
subplot(1,2,2);<br />
plot(pwm.time, pwm.motor_freq)<br />
hold on<br />
plot(rpm.time, rpm.motor_freq)<br />
plot(pwm.time(data_range), sim_rpm_freq)<br />
grid on<br />
xlim([0, max(pwm.time)]);<br />
ylim([0, 1.2*max(rpm.motor_freq)]);<br />
plot([t_min t_min], ylim, 'black');<br />
plot([t_max t_max], ylim, 'black');<br />
title(sprintf("PWM signal and response for \\tau = %.2f s", tau));<br />
xlabel("Time (s)");<br />
ylabel("Duty time (ms) / motor speed (rpm)");<br />
legend(["PWM signal", "RPM measurement", ...<br />
"Response approximation"])<br />
end<br />
<br />
end<br />
</source></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=File:Actuator_Dynamics_Measurements_Schematic.svg&diff=25210
File:Actuator Dynamics Measurements Schematic.svg
2020-10-18T20:46:18Z
<p>Hjkarssies: Hjkarssies uploaded a new version of File:Actuator Dynamics Measurements Schematic.svg</p>
<hr />
<div></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=File:Actuator_Dynamics_Measurements_Schematic.svg&diff=25209
File:Actuator Dynamics Measurements Schematic.svg
2020-10-18T20:13:17Z
<p>Hjkarssies: Hjkarssies uploaded a new version of File:Actuator Dynamics Measurements Schematic.svg</p>
<hr />
<div></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=File:Actuator_Dynamics_Measurements_Schematic.svg&diff=25208
File:Actuator Dynamics Measurements Schematic.svg
2020-10-18T20:10:16Z
<p>Hjkarssies: Hjkarssies uploaded a new version of File:Actuator Dynamics Measurements Schematic.svg</p>
<hr />
<div></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=File:Actuator_Dynamics_Measurements_Schematic.svg&diff=25207
File:Actuator Dynamics Measurements Schematic.svg
2020-10-18T20:09:28Z
<p>Hjkarssies: </p>
<hr />
<div></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=Actuator_Dynamics&diff=25206
Actuator Dynamics
2020-10-18T17:13:01Z
<p>Hjkarssies: </p>
<hr />
<div>[[File:Rpm_sensor.jpg|An RPM sensor|thumb]]<br />
[[File:Saleae-logic.jpg|A Saleae logic analyzer|thumb]]<br />
[[File:Saleae-analog-example.png|A screenshot from the Saleae logic analyser software|thumb]]<br />
[[File:Actuator_dynamics_matlab.svg|An example of estimated actuator dynamics in Matlab|thumb]]<br />
<br />
Some stabilization modules require an estimation of the current actuator states. Examples are the [[Working with INDI|INDI]] and INCA modules. When actuator feedback is not available, this estimation can be mathematical model of the actuators. For brushless motors, a first order transfer function often suffices as actuator model. This guide describes how to accurately determine a first order approximation for a motor with a propeller. Note that the actuator dynamics vary for each configuration of different motors, propellers and ESCs (Electronic Speed Controllers). Used for this example are:<br />
* [https://www.cobramotorsusa.com/motors-2814-16.html Cobra C-2814/16 brushless motor, KV=1050]<br />
* [https://www.banggood.com/1-Piece-HQProp-7040-7x4-CW-or-CCW-Propeller-Slow-Flyer-2-Blade-Green-Nylon-or-Glass-Fiber-Composite-For-RC-Airplane-Multirotors-p-1550754.html HQProp 7040 7x4 Propeller Slow Flyer 2 Blade Nylon/Glass Fiber Composite]<br />
* [http://kiss.flyduino.net/kiss-esc-24a-re/ KISS ESC 24A re – 32bit Electronic Speed Controller]<br><br><br />
<br />
==Actuator model==<br />
With this method, actuators are approximated by a first order model in de following form:<br><br><br />
<math>H_{act}=\frac{K}{\tau s+1}</math><br><br><br />
where <math>H_{act}</math> is the actuator transfer function in the Laplace domain, <math>K</math> is the system's gain, and <math>\tau</math> is the system's time constant.<br><br><br />
<br />
When using such a model on discrete time systems, an estimation of an actuator's current position (<math>\delta_{est}</math>) can be calculated as a function of the current actuator command (<math>\delta</math>) and the previous estimation of its position (<math>\delta_{prev}</math>):<br><br><br />
<math>\delta_{est}=\delta_{prev}+\alpha(\delta-\delta_{prev})</math><br><br />
where <math>\alpha=1-e^{-\tau\Delta t}</math><br><br />
where <math>\Delta t</math> is the duration of one time step<br><br><br />
<br />
To make a good approximation of actuator dynamics, one needs to determine appropriate values for <math>\tau</math> and <math>\alpha</math><br><br><br />
<br />
==Measuring Actuator Response==<br />
The first step in approximating actuator dynamics is to measure actuator response.<br />
<br />
<source lang="Matlab"><br />
%{<br />
This script is written for estimation of first order actuator dynamics for<br />
motors with rotors. It uses measurements of both a PWM input signal and the<br />
actuator's RPM response, and tries to fit a first order model to these data<br />
of de form:<br />
<br />
K <br />
-----------<br />
tau * s + 1<br />
<br />
This form can be used to estimate a actuator's state as a function of its<br />
previous state and its current input as follows:<br />
<br />
u_est = u_prev + alpha * (u - u_prev)<br />
where alpha = 1 - e^(tau * dt)<br />
<br />
© H.J. Karssies - Delft University of Technology<br />
%}<br />
<br />
% Close figures and clear workspace<br />
clc<br />
clear all<br />
close all<br />
<br />
% Settings<br />
data_file = 'motor_data_04.mat'; % Data file<br />
pwm.channel_no = 0; % Channel number of PWM input data<br />
pwm.period_min = 0.0005; % Minimum PWM period used for filtering<br />
rpm.channel_no = 1; % Channel number of RPM measurement data<br />
rpm.period_min = 0.00005; % Minimum RPM period used for filtering<br />
t_min = 5; % Starting time of relevant data (s)<br />
t_max = 10.4; % End time of relevant data (s)<br />
taus = 5:5:100; % List of values for tau to try (s)<br />
main_freq = 512; % Paparazzi main frequency<br />
<br />
% Load data<br />
raw_data = load(data_file);<br />
pwm = load_channel(pwm, raw_data);<br />
rpm = load_channel(rpm, raw_data);<br />
pwm.freq = length(pwm.time)/(pwm.time(end)-pwm.time(1));<br />
<br />
% Calculate motor frequencies<br />
pwm.motor_freq = pwm.period*1000000;<br />
rpm.motor_freq = 1./rpm.period;<br />
<br />
% Check performance for different values of tau<br />
for n = 1:length(taus)<br />
errors(n) = calc_performance(taus(n), pwm, rpm, t_min, t_max, false);<br />
end<br />
<br />
% Interpolate results and find minimum<br />
taus_int = taus(1):0.01:taus(end);<br />
errors_int = spline(taus, errors, taus_int);<br />
error_min = min(errors_int);<br />
tau_min = taus_int(find(errors_int == error_min));<br />
alpha = 1-exp(-tau_min/main_freq);<br />
<br />
% Plot result<br />
figure;<br />
set(gcf, 'Position', [100, 100, 1080, 540]);<br />
subplot(1,2,1);<br />
plot(taus_int, errors_int);<br />
hold on<br />
scatter(tau_min, error_min, 'black');<br />
grid on<br />
title("\tau estimation");<br />
xlabel("\tau (s)");<br />
ylabel("Error (rpm²)");<br />
text(tau_min, error_min+0.1*diff(ylim), ...<br />
[sprintf("\\tau = %.2f s", tau_min); ...<br />
sprintf("\\alpha = %.5f (%d Hz)", [alpha, main_freq])]);<br />
<br />
% Perform calculations with final tau and plot figures<br />
calc_performance(tau_min, pwm, rpm, t_min, t_max, true);<br />
<br />
% Load channel data and process it<br />
function channel = load_channel(channel, data)<br />
<br />
% Process channel<br />
channel_data = eval('data.digital_channel_' + string(channel.channel_no));<br />
channel.bit_1 = data.digital_channel_initial_bitstates(channel.channel_no+1);<br />
channel.period = channel_data(2-channel.bit_1:2:end)./data.digital_sample_rate_hz;<br />
channel.time = cumsum(channel_data./data.digital_sample_rate_hz);<br />
channel.time = channel.time(2-channel.bit_1:2:end);<br />
<br />
% Filter out invalid data<br />
channel.time = channel.time(channel.period>channel.period_min);<br />
channel.period = channel.period(channel.period>channel.period_min);<br />
<br />
end<br />
<br />
% Calculate performance with certain tau<br />
function error = calc_performance(tau, pwm, rpm, t_min, t_max, plot_fig)<br />
<br />
% Simulate RPM response<br />
motor = tf(tau,[1 tau]);<br />
<br />
% Run simulation at rounded PWM sample times<br />
bias = mean(pwm.motor_freq(1:10));<br />
dt = 1/round(pwm.freq);<br />
t = [0:length(pwm.time)-1]*dt;<br />
motor = tf(tau,[1 tau]);<br />
sim_pwm_freq = lsim(motor, pwm.motor_freq-bias, t) + bias;<br />
<br />
% Get the RPM requency values at the rounden PWM sample times<br />
int_rpm_freq = interp1(rpm.time,rpm.motor_freq,t,[],1)';<br />
int_rpm_freq(isinf(int_rpm_freq)) = 1;<br />
<br />
% Fit a simple model y = ax + b<br />
data_range = and(t>=t_min, t<=t_max);<br />
inputs = [ones(size(sim_pwm_freq(data_range))) ...<br />
sim_pwm_freq(data_range)];<br />
gain = inputs \ int_rpm_freq(data_range);<br />
sim_rpm_freq = inputs*gain;<br />
<br />
% Calculate squared error<br />
error = sum((int_rpm_freq(data_range)-sim_rpm_freq).^2);<br />
<br />
% Plot figure<br />
if plot_fig<br />
subplot(1,2,2);<br />
plot(pwm.time, pwm.motor_freq)<br />
hold on<br />
plot(rpm.time, rpm.motor_freq)<br />
plot(pwm.time(data_range), sim_rpm_freq)<br />
grid on<br />
xlim([0, max(pwm.time)]);<br />
ylim([0, 1.2*max(rpm.motor_freq)]);<br />
plot([t_min t_min], ylim, 'black');<br />
plot([t_max t_max], ylim, 'black');<br />
title(sprintf("PWM signal and response for \\tau = %.2f s", tau));<br />
xlabel("Time (s)");<br />
ylabel("Duty time (ms) / motor speed (rpm)");<br />
legend(["PWM signal", "RPM measurement", ...<br />
"Response approximation"])<br />
end<br />
<br />
end<br />
</source></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=Actuator_Dynamics&diff=25173
Actuator Dynamics
2020-09-23T19:33:06Z
<p>Hjkarssies: </p>
<hr />
<div>[[File:Rpm_sensor.jpg|An RPM sensor|thumb]]<br />
[[File:Saleae-logic.jpg|A Saleae logic analyzer|thumb]]<br />
[[File:Saleae-analog-example.png|A screenshot from the Saleae logic analyser software|thumb]]<br />
[[File:Actuator_dynamics_matlab.svg|An example of estimated actuator dynamics in Matlab|thumb]]<br />
<br />
==Work in progess==<br />
<br />
Some stabilization modules require an estimation of the current actuator states. Examples are the INDI and INCA modules. When actuator feedback is not available, this estimation can be mathematical model of the actuators. This guide describes how to accurately determine a first order approximation for a motor with a propeller.<br />
<br />
<source lang="Matlab"><br />
%{<br />
This script is written for estimation of first order actuator dynamics for<br />
motors with rotors. It uses measurements of both a PWM input signal and the<br />
actuator's RPM response, and tries to fit a first order model to these data<br />
of de form:<br />
<br />
K <br />
-----------<br />
tau * s + 1<br />
<br />
This form can be used to estimate a actuator's state as a function of its<br />
previous state and its current input as follows:<br />
<br />
u_est = u_prev + alpha * (u - u_prev)<br />
where alpha = 1 - e^(tau * dt)<br />
<br />
© H.J. Karssies - Delft University of Technology<br />
%}<br />
<br />
% Close figures and clear workspace<br />
clc<br />
clear all<br />
close all<br />
<br />
% Settings<br />
data_file = 'motor_data_04.mat'; % Data file<br />
pwm.channel_no = 0; % Channel number of PWM input data<br />
pwm.period_min = 0.0005; % Minimum PWM period used for filtering<br />
rpm.channel_no = 1; % Channel number of RPM measurement data<br />
rpm.period_min = 0.00005; % Minimum RPM period used for filtering<br />
t_min = 5; % Starting time of relevant data (s)<br />
t_max = 10.4; % End time of relevant data (s)<br />
taus = 5:5:100; % List of values for tau to try (s)<br />
main_freq = 512; % Paparazzi main frequency<br />
<br />
% Load data<br />
raw_data = load(data_file);<br />
pwm = load_channel(pwm, raw_data);<br />
rpm = load_channel(rpm, raw_data);<br />
pwm.freq = length(pwm.time)/(pwm.time(end)-pwm.time(1));<br />
<br />
% Calculate motor frequencies<br />
pwm.motor_freq = pwm.period*1000000;<br />
rpm.motor_freq = 1./rpm.period;<br />
<br />
% Check performance for different values of tau<br />
for n = 1:length(taus)<br />
errors(n) = calc_performance(taus(n), pwm, rpm, t_min, t_max, false);<br />
end<br />
<br />
% Interpolate results and find minimum<br />
taus_int = taus(1):0.01:taus(end);<br />
errors_int = spline(taus, errors, taus_int);<br />
error_min = min(errors_int);<br />
tau_min = taus_int(find(errors_int == error_min));<br />
alpha = 1-exp(-tau_min/main_freq);<br />
<br />
% Plot result<br />
figure;<br />
set(gcf, 'Position', [100, 100, 1080, 540]);<br />
subplot(1,2,1);<br />
plot(taus_int, errors_int);<br />
hold on<br />
scatter(tau_min, error_min, 'black');<br />
grid on<br />
title("\tau estimation");<br />
xlabel("\tau (s)");<br />
ylabel("Error (rpm²)");<br />
text(tau_min, error_min+0.1*diff(ylim), ...<br />
[sprintf("\\tau = %.2f s", tau_min); ...<br />
sprintf("\\alpha = %.5f (%d Hz)", [alpha, main_freq])]);<br />
<br />
% Perform calculations with final tau and plot figures<br />
calc_performance(tau_min, pwm, rpm, t_min, t_max, true);<br />
<br />
% Load channel data and process it<br />
function channel = load_channel(channel, data)<br />
<br />
% Process channel<br />
channel_data = eval('data.digital_channel_' + string(channel.channel_no));<br />
channel.bit_1 = data.digital_channel_initial_bitstates(channel.channel_no+1);<br />
channel.period = channel_data(2-channel.bit_1:2:end)./data.digital_sample_rate_hz;<br />
channel.time = cumsum(channel_data./data.digital_sample_rate_hz);<br />
channel.time = channel.time(2-channel.bit_1:2:end);<br />
<br />
% Filter out invalid data<br />
channel.time = channel.time(channel.period>channel.period_min);<br />
channel.period = channel.period(channel.period>channel.period_min);<br />
<br />
end<br />
<br />
% Calculate performance with certain tau<br />
function error = calc_performance(tau, pwm, rpm, t_min, t_max, plot_fig)<br />
<br />
% Simulate RPM response<br />
motor = tf(tau,[1 tau]);<br />
<br />
% Run simulation at rounded PWM sample times<br />
bias = mean(pwm.motor_freq(1:10));<br />
dt = 1/round(pwm.freq);<br />
t = [0:length(pwm.time)-1]*dt;<br />
motor = tf(tau,[1 tau]);<br />
sim_pwm_freq = lsim(motor, pwm.motor_freq-bias, t) + bias;<br />
<br />
% Get the RPM requency values at the rounden PWM sample times<br />
int_rpm_freq = interp1(rpm.time,rpm.motor_freq,t,[],1)';<br />
int_rpm_freq(isinf(int_rpm_freq)) = 1;<br />
<br />
% Fit a simple model y = ax + b<br />
data_range = and(t>=t_min, t<=t_max);<br />
inputs = [ones(size(sim_pwm_freq(data_range))) ...<br />
sim_pwm_freq(data_range)];<br />
gain = inputs \ int_rpm_freq(data_range);<br />
sim_rpm_freq = inputs*gain;<br />
<br />
% Calculate squared error<br />
error = sum((int_rpm_freq(data_range)-sim_rpm_freq).^2);<br />
<br />
% Plot figure<br />
if plot_fig<br />
subplot(1,2,2);<br />
plot(pwm.time, pwm.motor_freq)<br />
hold on<br />
plot(rpm.time, rpm.motor_freq)<br />
plot(pwm.time(data_range), sim_rpm_freq)<br />
grid on<br />
xlim([0, max(pwm.time)]);<br />
ylim([0, 1.2*max(rpm.motor_freq)]);<br />
plot([t_min t_min], ylim, 'black');<br />
plot([t_max t_max], ylim, 'black');<br />
title(sprintf("PWM signal and response for \\tau = %.2f s", tau));<br />
xlabel("Time (s)");<br />
ylabel("Duty time (ms) / motor speed (rpm)");<br />
legend(["PWM signal", "RPM measurement", ...<br />
"Response approximation"])<br />
end<br />
<br />
end<br />
</source></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=File:Saleae-analog-example.png&diff=25172
File:Saleae-analog-example.png
2020-09-23T19:27:38Z
<p>Hjkarssies: </p>
<hr />
<div></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=File:Saleae-logic.jpg&diff=25171
File:Saleae-logic.jpg
2020-09-23T19:27:20Z
<p>Hjkarssies: </p>
<hr />
<div></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=Actuator_Dynamics&diff=25170
Actuator Dynamics
2020-09-23T19:23:46Z
<p>Hjkarssies: </p>
<hr />
<div>[[File:Actuator_dynamics_matlab.svg|An example of estimated actuator dynamics in Matlab|500px|thumb]]<br />
[[File:Rpm_sensor.jpg|An RPM sensor|thumb]]<br />
<br />
==Work in progess==<br />
<br />
Some stabilization modules require an estimation of the current actuator states. Examples are the INDI and INCA modules. When actuator feedback is not available, this estimation can be mathematical model of the actuators. This guide describes how to accurately determine a first order approximation for a motor with a propeller.<br />
<br />
<source lang="Matlab"><br />
%{<br />
This script is written for estimation of first order actuator dynamics for<br />
motors with rotors. It uses measurements of both a PWM input signal and the<br />
actuator's RPM response, and tries to fit a first order model to these data<br />
of de form:<br />
<br />
K <br />
-----------<br />
tau * s + 1<br />
<br />
This form can be used to estimate a actuator's state as a function of its<br />
previous state and its current input as follows:<br />
<br />
u_est = u_prev + alpha * (u - u_prev)<br />
where alpha = 1 - e^(tau * dt)<br />
<br />
© H.J. Karssies - Delft University of Technology<br />
%}<br />
<br />
% Close figures and clear workspace<br />
clc<br />
clear all<br />
close all<br />
<br />
% Settings<br />
data_file = 'motor_data_04.mat'; % Data file<br />
pwm.channel_no = 0; % Channel number of PWM input data<br />
pwm.period_min = 0.0005; % Minimum PWM period used for filtering<br />
rpm.channel_no = 1; % Channel number of RPM measurement data<br />
rpm.period_min = 0.00005; % Minimum RPM period used for filtering<br />
t_min = 5; % Starting time of relevant data (s)<br />
t_max = 10.4; % End time of relevant data (s)<br />
taus = 5:5:100; % List of values for tau to try (s)<br />
main_freq = 512; % Paparazzi main frequency<br />
<br />
% Load data<br />
raw_data = load(data_file);<br />
pwm = load_channel(pwm, raw_data);<br />
rpm = load_channel(rpm, raw_data);<br />
pwm.freq = length(pwm.time)/(pwm.time(end)-pwm.time(1));<br />
<br />
% Calculate motor frequencies<br />
pwm.motor_freq = pwm.period*1000000;<br />
rpm.motor_freq = 1./rpm.period;<br />
<br />
% Check performance for different values of tau<br />
for n = 1:length(taus)<br />
errors(n) = calc_performance(taus(n), pwm, rpm, t_min, t_max, false);<br />
end<br />
<br />
% Interpolate results and find minimum<br />
taus_int = taus(1):0.01:taus(end);<br />
errors_int = spline(taus, errors, taus_int);<br />
error_min = min(errors_int);<br />
tau_min = taus_int(find(errors_int == error_min));<br />
alpha = 1-exp(-tau_min/main_freq);<br />
<br />
% Plot result<br />
figure;<br />
set(gcf, 'Position', [100, 100, 1080, 540]);<br />
subplot(1,2,1);<br />
plot(taus_int, errors_int);<br />
hold on<br />
scatter(tau_min, error_min, 'black');<br />
grid on<br />
title("\tau estimation");<br />
xlabel("\tau (s)");<br />
ylabel("Error (rpm²)");<br />
text(tau_min, error_min+0.1*diff(ylim), ...<br />
[sprintf("\\tau = %.2f s", tau_min); ...<br />
sprintf("\\alpha = %.5f (%d Hz)", [alpha, main_freq])]);<br />
<br />
% Perform calculations with final tau and plot figures<br />
calc_performance(tau_min, pwm, rpm, t_min, t_max, true);<br />
<br />
% Load channel data and process it<br />
function channel = load_channel(channel, data)<br />
<br />
% Process channel<br />
channel_data = eval('data.digital_channel_' + string(channel.channel_no));<br />
channel.bit_1 = data.digital_channel_initial_bitstates(channel.channel_no+1);<br />
channel.period = channel_data(2-channel.bit_1:2:end)./data.digital_sample_rate_hz;<br />
channel.time = cumsum(channel_data./data.digital_sample_rate_hz);<br />
channel.time = channel.time(2-channel.bit_1:2:end);<br />
<br />
% Filter out invalid data<br />
channel.time = channel.time(channel.period>channel.period_min);<br />
channel.period = channel.period(channel.period>channel.period_min);<br />
<br />
end<br />
<br />
% Calculate performance with certain tau<br />
function error = calc_performance(tau, pwm, rpm, t_min, t_max, plot_fig)<br />
<br />
% Simulate RPM response<br />
motor = tf(tau,[1 tau]);<br />
<br />
% Run simulation at rounded PWM sample times<br />
bias = mean(pwm.motor_freq(1:10));<br />
dt = 1/round(pwm.freq);<br />
t = [0:length(pwm.time)-1]*dt;<br />
motor = tf(tau,[1 tau]);<br />
sim_pwm_freq = lsim(motor, pwm.motor_freq-bias, t) + bias;<br />
<br />
% Get the RPM requency values at the rounden PWM sample times<br />
int_rpm_freq = interp1(rpm.time,rpm.motor_freq,t,[],1)';<br />
int_rpm_freq(isinf(int_rpm_freq)) = 1;<br />
<br />
% Fit a simple model y = ax + b<br />
data_range = and(t>=t_min, t<=t_max);<br />
inputs = [ones(size(sim_pwm_freq(data_range))) ...<br />
sim_pwm_freq(data_range)];<br />
gain = inputs \ int_rpm_freq(data_range);<br />
sim_rpm_freq = inputs*gain;<br />
<br />
% Calculate squared error<br />
error = sum((int_rpm_freq(data_range)-sim_rpm_freq).^2);<br />
<br />
% Plot figure<br />
if plot_fig<br />
subplot(1,2,2);<br />
plot(pwm.time, pwm.motor_freq)<br />
hold on<br />
plot(rpm.time, rpm.motor_freq)<br />
plot(pwm.time(data_range), sim_rpm_freq)<br />
grid on<br />
xlim([0, max(pwm.time)]);<br />
ylim([0, 1.2*max(rpm.motor_freq)]);<br />
plot([t_min t_min], ylim, 'black');<br />
plot([t_max t_max], ylim, 'black');<br />
title(sprintf("PWM signal and response for \\tau = %.2f s", tau));<br />
xlabel("Time (s)");<br />
ylabel("Duty time (ms) / motor speed (rpm)");<br />
legend(["PWM signal", "RPM measurement", ...<br />
"Response approximation"])<br />
end<br />
<br />
end<br />
</source></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=File:Rpm_sensor.jpg&diff=25169
File:Rpm sensor.jpg
2020-09-23T19:21:24Z
<p>Hjkarssies: </p>
<hr />
<div></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=Actuator_Dynamics&diff=25168
Actuator Dynamics
2020-09-23T19:19:08Z
<p>Hjkarssies: </p>
<hr />
<div>Work in progess<br />
<br />
Some stabilization modules require an estimation of the current actuator states. Examples are the INDI and INCA modules. When actuator feedback is not available, this estimation can be mathematical model of the actuators. This guide describes how to accurately determine a first order approximation for a motor with a propeller.<br />
<br />
[[File:Actuator_dynamics_matlab.svg|An example of estimated actuator dynamics in Matlab|frame|right|300px]]<br />
<br />
<source lang="Matlab"><br />
%{<br />
This script is written for estimation of first order actuator dynamics for<br />
motors with rotors. It uses measurements of both a PWM input signal and the<br />
actuator's RPM response, and tries to fit a first order model to these data<br />
of de form:<br />
<br />
K <br />
-----------<br />
tau * s + 1<br />
<br />
This form can be used to estimate a actuator's state as a function of its<br />
previous state and its current input as follows:<br />
<br />
u_est = u_prev + alpha * (u - u_prev)<br />
where alpha = 1 - e^(tau * dt)<br />
<br />
© H.J. Karssies - Delft University of Technology<br />
%}<br />
<br />
% Close figures and clear workspace<br />
clc<br />
clear all<br />
close all<br />
<br />
% Settings<br />
data_file = 'motor_data_04.mat'; % Data file<br />
pwm.channel_no = 0; % Channel number of PWM input data<br />
pwm.period_min = 0.0005; % Minimum PWM period used for filtering<br />
rpm.channel_no = 1; % Channel number of RPM measurement data<br />
rpm.period_min = 0.00005; % Minimum RPM period used for filtering<br />
t_min = 5; % Starting time of relevant data (s)<br />
t_max = 10.4; % End time of relevant data (s)<br />
taus = 5:5:100; % List of values for tau to try (s)<br />
main_freq = 512; % Paparazzi main frequency<br />
<br />
% Load data<br />
raw_data = load(data_file);<br />
pwm = load_channel(pwm, raw_data);<br />
rpm = load_channel(rpm, raw_data);<br />
pwm.freq = length(pwm.time)/(pwm.time(end)-pwm.time(1));<br />
<br />
% Calculate motor frequencies<br />
pwm.motor_freq = pwm.period*1000000;<br />
rpm.motor_freq = 1./rpm.period;<br />
<br />
% Check performance for different values of tau<br />
for n = 1:length(taus)<br />
errors(n) = calc_performance(taus(n), pwm, rpm, t_min, t_max, false);<br />
end<br />
<br />
% Interpolate results and find minimum<br />
taus_int = taus(1):0.01:taus(end);<br />
errors_int = spline(taus, errors, taus_int);<br />
error_min = min(errors_int);<br />
tau_min = taus_int(find(errors_int == error_min));<br />
alpha = 1-exp(-tau_min/main_freq);<br />
<br />
% Plot result<br />
figure;<br />
set(gcf, 'Position', [100, 100, 1080, 540]);<br />
subplot(1,2,1);<br />
plot(taus_int, errors_int);<br />
hold on<br />
scatter(tau_min, error_min, 'black');<br />
grid on<br />
title("\tau estimation");<br />
xlabel("\tau (s)");<br />
ylabel("Error (rpm²)");<br />
text(tau_min, error_min+0.1*diff(ylim), ...<br />
[sprintf("\\tau = %.2f s", tau_min); ...<br />
sprintf("\\alpha = %.5f (%d Hz)", [alpha, main_freq])]);<br />
<br />
% Perform calculations with final tau and plot figures<br />
calc_performance(tau_min, pwm, rpm, t_min, t_max, true);<br />
<br />
% Load channel data and process it<br />
function channel = load_channel(channel, data)<br />
<br />
% Process channel<br />
channel_data = eval('data.digital_channel_' + string(channel.channel_no));<br />
channel.bit_1 = data.digital_channel_initial_bitstates(channel.channel_no+1);<br />
channel.period = channel_data(2-channel.bit_1:2:end)./data.digital_sample_rate_hz;<br />
channel.time = cumsum(channel_data./data.digital_sample_rate_hz);<br />
channel.time = channel.time(2-channel.bit_1:2:end);<br />
<br />
% Filter out invalid data<br />
channel.time = channel.time(channel.period>channel.period_min);<br />
channel.period = channel.period(channel.period>channel.period_min);<br />
<br />
end<br />
<br />
% Calculate performance with certain tau<br />
function error = calc_performance(tau, pwm, rpm, t_min, t_max, plot_fig)<br />
<br />
% Simulate RPM response<br />
motor = tf(tau,[1 tau]);<br />
<br />
% Run simulation at rounded PWM sample times<br />
bias = mean(pwm.motor_freq(1:10));<br />
dt = 1/round(pwm.freq);<br />
t = [0:length(pwm.time)-1]*dt;<br />
motor = tf(tau,[1 tau]);<br />
sim_pwm_freq = lsim(motor, pwm.motor_freq-bias, t) + bias;<br />
<br />
% Get the RPM requency values at the rounden PWM sample times<br />
int_rpm_freq = interp1(rpm.time,rpm.motor_freq,t,[],1)';<br />
int_rpm_freq(isinf(int_rpm_freq)) = 1;<br />
<br />
% Fit a simple model y = ax + b<br />
data_range = and(t>=t_min, t<=t_max);<br />
inputs = [ones(size(sim_pwm_freq(data_range))) ...<br />
sim_pwm_freq(data_range)];<br />
gain = inputs \ int_rpm_freq(data_range);<br />
sim_rpm_freq = inputs*gain;<br />
<br />
% Calculate squared error<br />
error = sum((int_rpm_freq(data_range)-sim_rpm_freq).^2);<br />
<br />
% Plot figure<br />
if plot_fig<br />
subplot(1,2,2);<br />
plot(pwm.time, pwm.motor_freq)<br />
hold on<br />
plot(rpm.time, rpm.motor_freq)<br />
plot(pwm.time(data_range), sim_rpm_freq)<br />
grid on<br />
xlim([0, max(pwm.time)]);<br />
ylim([0, 1.2*max(rpm.motor_freq)]);<br />
plot([t_min t_min], ylim, 'black');<br />
plot([t_max t_max], ylim, 'black');<br />
title(sprintf("PWM signal and response for \\tau = %.2f s", tau));<br />
xlabel("Time (s)");<br />
ylabel("Duty time (ms) / motor speed (rpm)");<br />
legend(["PWM signal", "RPM measurement", ...<br />
"Response approximation"])<br />
end<br />
<br />
end<br />
</source></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=Actuator_Dynamics&diff=25167
Actuator Dynamics
2020-09-23T19:18:21Z
<p>Hjkarssies: </p>
<hr />
<div>Work in progess<br />
<br />
Some stabilization modules require an estimation of the current actuator states. Examples are the INDI and INCA modules. When actuator feedback is not available, this estimation can be mathematical model of the actuators. This guide describes how to accurately determine a first order approximation for a motor with a propeller.<br />
<br />
[[File:Actuator_dynamics_matlab.svg|An example of estimated actuator dynamics in Matlab]]<br />
<br />
<source lang="Matlab"><br />
%{<br />
This script is written for estimation of first order actuator dynamics for<br />
motors with rotors. It uses measurements of both a PWM input signal and the<br />
actuator's RPM response, and tries to fit a first order model to these data<br />
of de form:<br />
<br />
K <br />
-----------<br />
tau * s + 1<br />
<br />
This form can be used to estimate a actuator's state as a function of its<br />
previous state and its current input as follows:<br />
<br />
u_est = u_prev + alpha * (u - u_prev)<br />
where alpha = 1 - e^(tau * dt)<br />
<br />
© H.J. Karssies - Delft University of Technology<br />
%}<br />
<br />
% Close figures and clear workspace<br />
clc<br />
clear all<br />
close all<br />
<br />
% Settings<br />
data_file = 'motor_data_04.mat'; % Data file<br />
pwm.channel_no = 0; % Channel number of PWM input data<br />
pwm.period_min = 0.0005; % Minimum PWM period used for filtering<br />
rpm.channel_no = 1; % Channel number of RPM measurement data<br />
rpm.period_min = 0.00005; % Minimum RPM period used for filtering<br />
t_min = 5; % Starting time of relevant data (s)<br />
t_max = 10.4; % End time of relevant data (s)<br />
taus = 5:5:100; % List of values for tau to try (s)<br />
main_freq = 512; % Paparazzi main frequency<br />
<br />
% Load data<br />
raw_data = load(data_file);<br />
pwm = load_channel(pwm, raw_data);<br />
rpm = load_channel(rpm, raw_data);<br />
pwm.freq = length(pwm.time)/(pwm.time(end)-pwm.time(1));<br />
<br />
% Calculate motor frequencies<br />
pwm.motor_freq = pwm.period*1000000;<br />
rpm.motor_freq = 1./rpm.period;<br />
<br />
% Check performance for different values of tau<br />
for n = 1:length(taus)<br />
errors(n) = calc_performance(taus(n), pwm, rpm, t_min, t_max, false);<br />
end<br />
<br />
% Interpolate results and find minimum<br />
taus_int = taus(1):0.01:taus(end);<br />
errors_int = spline(taus, errors, taus_int);<br />
error_min = min(errors_int);<br />
tau_min = taus_int(find(errors_int == error_min));<br />
alpha = 1-exp(-tau_min/main_freq);<br />
<br />
% Plot result<br />
figure;<br />
set(gcf, 'Position', [100, 100, 1080, 540]);<br />
subplot(1,2,1);<br />
plot(taus_int, errors_int);<br />
hold on<br />
scatter(tau_min, error_min, 'black');<br />
grid on<br />
title("\tau estimation");<br />
xlabel("\tau (s)");<br />
ylabel("Error (rpm²)");<br />
text(tau_min, error_min+0.1*diff(ylim), ...<br />
[sprintf("\\tau = %.2f s", tau_min); ...<br />
sprintf("\\alpha = %.5f (%d Hz)", [alpha, main_freq])]);<br />
<br />
% Perform calculations with final tau and plot figures<br />
calc_performance(tau_min, pwm, rpm, t_min, t_max, true);<br />
<br />
% Load channel data and process it<br />
function channel = load_channel(channel, data)<br />
<br />
% Process channel<br />
channel_data = eval('data.digital_channel_' + string(channel.channel_no));<br />
channel.bit_1 = data.digital_channel_initial_bitstates(channel.channel_no+1);<br />
channel.period = channel_data(2-channel.bit_1:2:end)./data.digital_sample_rate_hz;<br />
channel.time = cumsum(channel_data./data.digital_sample_rate_hz);<br />
channel.time = channel.time(2-channel.bit_1:2:end);<br />
<br />
% Filter out invalid data<br />
channel.time = channel.time(channel.period>channel.period_min);<br />
channel.period = channel.period(channel.period>channel.period_min);<br />
<br />
end<br />
<br />
% Calculate performance with certain tau<br />
function error = calc_performance(tau, pwm, rpm, t_min, t_max, plot_fig)<br />
<br />
% Simulate RPM response<br />
motor = tf(tau,[1 tau]);<br />
<br />
% Run simulation at rounded PWM sample times<br />
bias = mean(pwm.motor_freq(1:10));<br />
dt = 1/round(pwm.freq);<br />
t = [0:length(pwm.time)-1]*dt;<br />
motor = tf(tau,[1 tau]);<br />
sim_pwm_freq = lsim(motor, pwm.motor_freq-bias, t) + bias;<br />
<br />
% Get the RPM requency values at the rounden PWM sample times<br />
int_rpm_freq = interp1(rpm.time,rpm.motor_freq,t,[],1)';<br />
int_rpm_freq(isinf(int_rpm_freq)) = 1;<br />
<br />
% Fit a simple model y = ax + b<br />
data_range = and(t>=t_min, t<=t_max);<br />
inputs = [ones(size(sim_pwm_freq(data_range))) ...<br />
sim_pwm_freq(data_range)];<br />
gain = inputs \ int_rpm_freq(data_range);<br />
sim_rpm_freq = inputs*gain;<br />
<br />
% Calculate squared error<br />
error = sum((int_rpm_freq(data_range)-sim_rpm_freq).^2);<br />
<br />
% Plot figure<br />
if plot_fig<br />
subplot(1,2,2);<br />
plot(pwm.time, pwm.motor_freq)<br />
hold on<br />
plot(rpm.time, rpm.motor_freq)<br />
plot(pwm.time(data_range), sim_rpm_freq)<br />
grid on<br />
xlim([0, max(pwm.time)]);<br />
ylim([0, 1.2*max(rpm.motor_freq)]);<br />
plot([t_min t_min], ylim, 'black');<br />
plot([t_max t_max], ylim, 'black');<br />
title(sprintf("PWM signal and response for \\tau = %.2f s", tau));<br />
xlabel("Time (s)");<br />
ylabel("Duty time (ms) / motor speed (rpm)");<br />
legend(["PWM signal", "RPM measurement", ...<br />
"Response approximation"])<br />
end<br />
<br />
end<br />
</source></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=File:Actuator_dynamics_matlab.svg&diff=25166
File:Actuator dynamics matlab.svg
2020-09-23T19:17:20Z
<p>Hjkarssies: </p>
<hr />
<div></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=Actuator_Dynamics&diff=25165
Actuator Dynamics
2020-09-23T19:05:28Z
<p>Hjkarssies: Added code</p>
<hr />
<div>Work in progess<br />
<br />
Some stabilization modules require an estimation of the current actuator states. Examples are the INDI and INCA modules. When actuator feedback is not available, this estimation can be mathematical model of the actuators. This guide describes how to accurately determine a first order approximation for a motor with a propeller.<br />
<br />
<source lang="Matlab"><br />
%{<br />
This script is written for estimation of first order actuator dynamics for<br />
motors with rotors. It uses measurements of both a PWM input signal and the<br />
actuator's RPM response, and tries to fit a first order model to these data<br />
of de form:<br />
<br />
K <br />
-----------<br />
tau * s + 1<br />
<br />
This form can be used to estimate a actuator's state as a function of its<br />
previous state and its current input as follows:<br />
<br />
u_est = u_prev + alpha * (u - u_prev)<br />
where alpha = 1 - e^(tau * dt)<br />
<br />
© H.J. Karssies - Delft University of Technology<br />
%}<br />
<br />
% Close figures and clear workspace<br />
clc<br />
clear all<br />
close all<br />
<br />
% Settings<br />
data_file = 'motor_data_04.mat'; % Data file<br />
pwm.channel_no = 0; % Channel number of PWM input data<br />
pwm.period_min = 0.0005; % Minimum PWM period used for filtering<br />
rpm.channel_no = 1; % Channel number of RPM measurement data<br />
rpm.period_min = 0.00005; % Minimum RPM period used for filtering<br />
t_min = 5; % Starting time of relevant data (s)<br />
t_max = 10.4; % End time of relevant data (s)<br />
taus = 5:5:100; % List of values for tau to try (s)<br />
main_freq = 512; % Paparazzi main frequency<br />
<br />
% Load data<br />
raw_data = load(data_file);<br />
pwm = load_channel(pwm, raw_data);<br />
rpm = load_channel(rpm, raw_data);<br />
pwm.freq = length(pwm.time)/(pwm.time(end)-pwm.time(1));<br />
<br />
% Calculate motor frequencies<br />
pwm.motor_freq = pwm.period*1000000;<br />
rpm.motor_freq = 1./rpm.period;<br />
<br />
% Check performance for different values of tau<br />
for n = 1:length(taus)<br />
errors(n) = calc_performance(taus(n), pwm, rpm, t_min, t_max, false);<br />
end<br />
<br />
% Interpolate results and find minimum<br />
taus_int = taus(1):0.01:taus(end);<br />
errors_int = spline(taus, errors, taus_int);<br />
error_min = min(errors_int);<br />
tau_min = taus_int(find(errors_int == error_min));<br />
alpha = 1-exp(-tau_min/main_freq);<br />
<br />
% Plot result<br />
figure;<br />
set(gcf, 'Position', [100, 100, 1080, 540]);<br />
subplot(1,2,1);<br />
plot(taus_int, errors_int);<br />
hold on<br />
scatter(tau_min, error_min, 'black');<br />
grid on<br />
title("\tau estimation");<br />
xlabel("\tau (s)");<br />
ylabel("Error (rpm²)");<br />
text(tau_min, error_min+0.1*diff(ylim), ...<br />
[sprintf("\\tau = %.2f s", tau_min); ...<br />
sprintf("\\alpha = %.5f (%d Hz)", [alpha, main_freq])]);<br />
<br />
% Perform calculations with final tau and plot figures<br />
calc_performance(tau_min, pwm, rpm, t_min, t_max, true);<br />
<br />
% Load channel data and process it<br />
function channel = load_channel(channel, data)<br />
<br />
% Process channel<br />
channel_data = eval('data.digital_channel_' + string(channel.channel_no));<br />
channel.bit_1 = data.digital_channel_initial_bitstates(channel.channel_no+1);<br />
channel.period = channel_data(2-channel.bit_1:2:end)./data.digital_sample_rate_hz;<br />
channel.time = cumsum(channel_data./data.digital_sample_rate_hz);<br />
channel.time = channel.time(2-channel.bit_1:2:end);<br />
<br />
% Filter out invalid data<br />
channel.time = channel.time(channel.period>channel.period_min);<br />
channel.period = channel.period(channel.period>channel.period_min);<br />
<br />
end<br />
<br />
% Calculate performance with certain tau<br />
function error = calc_performance(tau, pwm, rpm, t_min, t_max, plot_fig)<br />
<br />
% Simulate RPM response<br />
motor = tf(tau,[1 tau]);<br />
<br />
% Run simulation at rounded PWM sample times<br />
bias = mean(pwm.motor_freq(1:10));<br />
dt = 1/round(pwm.freq);<br />
t = [0:length(pwm.time)-1]*dt;<br />
motor = tf(tau,[1 tau]);<br />
sim_pwm_freq = lsim(motor, pwm.motor_freq-bias, t) + bias;<br />
<br />
% Get the RPM requency values at the rounden PWM sample times<br />
int_rpm_freq = interp1(rpm.time,rpm.motor_freq,t,[],1)';<br />
int_rpm_freq(isinf(int_rpm_freq)) = 1;<br />
<br />
% Fit a simple model y = ax + b<br />
data_range = and(t>=t_min, t<=t_max);<br />
inputs = [ones(size(sim_pwm_freq(data_range))) ...<br />
sim_pwm_freq(data_range)];<br />
gain = inputs \ int_rpm_freq(data_range);<br />
sim_rpm_freq = inputs*gain;<br />
<br />
% Calculate squared error<br />
error = sum((int_rpm_freq(data_range)-sim_rpm_freq).^2);<br />
<br />
% Plot figure<br />
if plot_fig<br />
subplot(1,2,2);<br />
plot(pwm.time, pwm.motor_freq)<br />
hold on<br />
plot(rpm.time, rpm.motor_freq)<br />
plot(pwm.time(data_range), sim_rpm_freq)<br />
grid on<br />
xlim([0, max(pwm.time)]);<br />
ylim([0, 1.2*max(rpm.motor_freq)]);<br />
plot([t_min t_min], ylim, 'black');<br />
plot([t_max t_max], ylim, 'black');<br />
title(sprintf("PWM signal and response for \\tau = %.2f s", tau));<br />
xlabel("Time (s)");<br />
ylabel("Duty time (ms) / motor speed (rpm)");<br />
legend(["PWM signal", "RPM measurement", ...<br />
"Response approximation"])<br />
end<br />
<br />
end<br />
</source></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=Actuator_Dynamics&diff=25164
Actuator Dynamics
2020-09-23T19:05:04Z
<p>Hjkarssies: </p>
<hr />
<div>Work in progess<br />
<br />
Some stabilization modules require an estimation of the current actuator states. Examples are the INDI and INCA modules. When actuator feedback is not available, this estimation can be mathematical model of the actuators. This guide describes how to accurately determine a first order approximation for a motor with a propeller.<br />
<br />
<source lang="Matlab"><br />
Test<br />
</source></div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=Actuator_Dynamics&diff=25163
Actuator Dynamics
2020-09-23T19:03:33Z
<p>Hjkarssies: </p>
<hr />
<div>Work in progess<br />
<br />
Some stabilization modules require an estimation of the current actuator states. Examples are the INDI and INCA modules. When actuator feedback is not available, this estimation can be mathematical model of the actuators. This guide describes how to accurately determine a first order approximation for a motor with a propeller.<br />
<br />
{{inline-code|test}}</div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=Actuator_Dynamics&diff=25162
Actuator Dynamics
2020-09-23T19:03:00Z
<p>Hjkarssies: </p>
<hr />
<div>Work in progess<br />
<br />
Some stabilization modules require an estimation of the current actuator states. Examples are the INDI and INCA modules. When actuator feedback is not available, this estimation can be mathematical model of the actuators. This guide describes how to accurately determine a first order approximation for a motor with a propeller.<br />
<br />
{{inline-code test}}</div>
Hjkarssies
http://wiki.paparazziuav.org/w/index.php?title=Actuator_Dynamics&diff=25150
Actuator Dynamics
2020-09-23T07:48:34Z
<p>Hjkarssies: Created page with "Work in progess Some stabilization modules require an estimation of the current actuator states. Examples are the INDI and INCA modules. When actuator feedback is not availab..."</p>
<hr />
<div>Work in progess<br />
<br />
Some stabilization modules require an estimation of the current actuator states. Examples are the INDI and INCA modules. When actuator feedback is not available, this estimation can be mathematical model of the actuators. This guide describes how to accurately determine a first order approximation for a motor with a propeller.</div>
Hjkarssies