Simple PID Controller with Simulink

I want to develop a simple PID controller for hovering, please see the attachment. The mass M (10 lbs) will start at h0 = 0 and the goal is to hover at h (6’ for example). The input to the plant is the throttle (from 0 to 100%, 0 to 20 lbs thrust) and the output is the altitude. So the desired altitude is the reference value and the equation of motion is simply f=ma. Could someone help me with the Simulink diagram? This not a homework problem, I'm helping to develop a VTOL propulsion system with undergrad students at USC and this would help us to jumpstart the project.

[Jeff Dillon]

Revisiting this question, if appropriate here. It appears that the solution assumes an instant response of the lift force. In reality, the lift is provided by a RC jet engine. The throttle response for a jet engine is linear from 0 to 5 seconds. That is, to go from zero thrust to full throttle takes 5 seconds. How would this affect this problem? Also, on this line, how was the 4.44482... value calculated? I want to confirm my understanding as the original mass and available lift force will be input parameters.

Lmax = 20*4.4482216153;

[Aquatris]

The 0-5 seconds can be accounted by a simple Rate Limiter block to prevent instantenous changes to the control signal, but rather giving it a maximum slope. This would mean your system is now less responsive, hence your controller performance would likely degrade and you might wanna retune your controller.

4.4482216153 is the constant used to convert lbs to Newton.

[Jeff Dillon]

@Aquatris Thanks! Let's assume that the dT is .02 second (step time in the model, real time). The change in thrust of the jet turbine is 14 pounds/second or .28 pounds in .02 second. What would be the rising slew rate for a Rate Limiter block? Also, I caught a major mistake in my plant model. I was assuming for some reason that the control input was throttle %. It's not at all, it's Newton force, I need to convert that into throttle % inside my hardware plant.

[Aquatris]

If you choose the sample time as inherited in the rate limiter block, you should be able to simply enter 14 as your slew rate.

My advice to you: To have a better understanding of these blocks, you might wanna create a simple model with step input, rate limiter and scope blocks. Play with the model, change the time step, change slew rate and see the effects. This way it would be more clear for you.

[Jeff Dillon]

@Aquatris Great suggestion. Also I'll probably need to convert to Newtons so 14 * 4.45. And this block would go just before the transfer function?

[Aquatris]

That kinda depends on the model of your system. For linear system with sufficiently fast response, the difference might be negligible. But there might be cases where the difference is significant, as when you put the block before the transfer function, the states of your system moves faster than the observed output (e.g., you apply step input to the system, the system output becomes 10 in 0.1 sec, the output takes 5 sec to reach that equilibrium due to rate limiter, you send 0 input 2 sec after the step input, at this point the rate limiter output value is around 4, and transfer function output is 10. within that 0.1 sec your transfer function takes to reach equlibrium, the observed output will still rise until your transfer function output is less than the observed output. Is this desired behaviour? If so then you should put it after transfer function, however if it is not, then you should put it before transfer function)

Answer 1

[Sam Chak]

Hi @Jeff Dillon

Not a VTOL person, but from the free-body diagram, Newton's 2nd Law can be applied with the basic assumptions of relatively low altitude compared to the Earth radius, and air resistance is neglected. The basic equation of motion is given by

where the weight of VTOL system, . Rearranging that

.

Since the accelation of the VTOL is the second-order derivative of the altitude with respect to time, , then

.

Taking the Laplace transform with zero initial conditions, you get

which can be rearranged to obtain the transfer function from plant inputs to plant output

Here, is the command input and is the disturbance.

If can be accurately estimated, then you can design

to counter the effect of the disturbance, and is the PID thing that you need to design

m    = 10*0.45359237;       % mass of VTOL
g    = -9.80665;            % gravitational acceleration
W    = m*g                  % weight of VTOL

W = -44.4822

Gp   = tf(1, [m 0 0])

Gp =
 
      1
  ---------
  4.536 s^2
 
Continuous-time transfer function.

Lmax = 20*4.4482216153;     % max lift force available
Umax = Lmax + W             % design limit for U(s) or output of Gu(s)

Umax = 44.4822

% Design of U(s)
Gc   = pidtune(Gp, 'PIDF', 0.95)

Gc =
 
             1            s    
  Kp + Ki * --- + Kd * --------
             s          Tf*s+1 

  with Kp = 0.76, Ki = 0.0338, Kd = 4.2, Tf = 0.193
 
Continuous-time PIDF controller in parallel form.

Gcl  = minreal(feedback(Gc*Gp, 1))

Gcl =
 
           4.961 s^2 + 0.8757 s + 0.03861
  ------------------------------------------------
  s^4 + 5.181 s^3 + 4.961 s^2 + 0.8757 s + 0.03861
 
Continuous-time transfer function.

Gu   = minreal(feedback(Gc, Gp))

Gu =
 
         22.5 s^4 + 3.972 s^3 + 0.1752 s^2
  ------------------------------------------------
  s^4 + 5.181 s^3 + 4.961 s^2 + 0.8757 s + 0.03861
 
Continuous-time transfer function.

t    = linspace(0, 50, 5001);
hd   = 6*0.3048;            % desired altitude
u    = hd*heaviside(t);
% Check if output of Gcl tracks hd
lsim(Gcl, u, t), grid on

% Check if output of Gu exceeds 44.4822
lsim(Gu, u, t/10), grid on

The Simulink model should look like this:

[Jeff Dillon]

This works perfect! Now if I want to test this on real hardware which outputs the altitude and accepts the lift force (via ROS2), where would I place these two blocks?

[Sam Chak]

Good to hear that, @Jeff Dillon. If you find the solution in MATLAB/Simulink is helpful, please consider accepting ✔ and voting 👍 the Answer. Thanks!

Regarding the ROS 2, I think you need to ask the hardware experts or MathWorks staff by opening a new question.

[Jeff Dillon]

Sorry for the confusion. I already have the ROS2 blocks and the hardware working, you can ignore that part. I was just wondering where in your model I would interface to the hardware. Where would I add the real-world input/outputs.

[Paul]

@Sam Chak

heavisde(0) = 1/2 (unless changed via sympref), which probably isn't what's desired in this example. I don't think it changes any of the conclusions, but thought I'd point it out because it does affect the responses at t = 0 and shortly after.

[Sam Chak]

Thanks @Paul for the highlight.

OP can replace these 3 lines

u    = hd*heaviside(t);
lsim(Gcl, u, t), grid on
lsim(Gu, u, t/10), grid on

with

opt  = stepDataOptions('StepAmplitude', hd);
step(Gcl, t, opt), grid on
step(Gu, t/10, opt), grid on

or replace only

u    = hd*heaviside(t);

with

Gr     = tf(hd, 1);
[u, t] = step(Gr, t);

@Jeff Dillon, regarding the interface, perhap these 3 articles are useful for you.

https://www.mathworks.com/discovery/hardware-in-the-loop-hil.html

https://www.mathworks.com/help/simscape/ug/what-is-hardware-in-the-loop-simulation.html

https://www.mathworks.com/help/slrealtime/gs/real-time-simulation-and-testing.html

[Jeff Dillon]

The insight here as been amazing, thanks! Yes, I've done real-time HIL in the past including this fun project https://aeropendulum.arizona.edu/ . In my case, can I just replace the Transfer Function with the actual Plant block?

[Sam Chak]

Hi @Jeff Dillon

Thanks for showing the AeroPendulum Profject and the applications of MATLAB. 👍 +1

Yes, you can replace it with actual Plant block. And you may need to tune the PID to get the satisfactory performance.

[Jeff Dillon]

I have a completed demo with both simulation and the plant, how do I embed quality images inline?

[Sam Chak]

Click the "photo" icon (in the INSERT tab), if you want to embed the image into this forum thread.

There are display settings for you to adjust the size of the displayed image.

[Jeff Dillon]

Here is my first draft of the implementation, it will be interesting to compare actual vs simulated. I'm using ROS2 and UDP to compare latency. If ROS2 proves fast enough, I'll use it throughout.

[Jeff Dillon]

@Paul @Sam Chak On this line, how was the 4.44482... value calculated? I want to confirm my understanding as the original mass and available lift force will be input parameters.

Lmax = 20*4.4482216153;

[Sam Chak]

Hi @Jeff Dillon,

Sorry for the late reply.

After checking the page, it is the unit conversion of the pound-force to Newton from the available 20 lbs thrust.

[Jeff Dillon]

Ah, got it! Thanks @Sam Chak One last question. The transfer function currently assumes an instant response time for the lifting force. In reality, the jet turbine takes 5 seconds to go from 0% to 100% throttle, linearly. How would I model this delay?

[Sam Chak]

Hi @Jeff Dillon

I'm no expert in delayed systems. Your interpretation of time-delay may be different from mine.

To approximately describe the transient period of 5 seconds (from 0 to almost 100%) in a mathematical model, you really need to find out whether the jet turbine (actuator) behaves linearly or nonlinearly. If it behaves like a first-order linear system, then the model in time-domain looks like the following.

Note that I didn't use the transfer function because I don't want you to get confused.

tspan  = linspace(0, 10, 10001);
x0     = 0;
[t, x] = ode45(@odefcn, tspan, x0);

% Solution Plot
plot(t, x, 'linewidth', 1.5), 
grid on, xlabel('t'), ylabel('x')

% check the output of jet turbine at 5 seconds.
idx = find(t == 5);     % find the location of t = 5
jet_output = x(idx)     % insert this location to the function almost like f(t)

jet_output = 0.9975

% 1st-order linear system
function xdot = odefcn(t, x)
    u    = 1;       % input reference
    Ts   = 5;       % transient time
    tau  = Ts/6;    % time constant parameter in the model
    xdot = (- x + u)/tau;
end

[Jeff Dillon]

@Sam Chak As mentioned, the response of the jet turbine is linear. At time = 0, if you command full throttle, it would take 5 seconds to get to full throttle, linearly. At 2.5 seconds, it would be at 1/2 throttle, etc. If the turbine is at full power, and you command it to zero throttle, it would take 5 sconds to throttle down, etc. What would the transfer function look like in the Simulink model? In your first graph, shouldn't it be a straight line? I guess I'm not sure what the y-axis represents. If it's the thrust force normalized between 0 and 1, it should indeed be a straight line.

[Sam Chak]

@Jeff Dillon,

I think that your description about the linear relationship of the input-output response of the throttle is actually a nonlinear system. For nonlinear systems, there is no transfer function.

To find out whether a system is linear or nonlinear, we look at the ODE, not the input-output response.

I'm not a jet turbine expert. Can you find out if the throttle in the jet turbine really behaves as linearly as you described?

[Jeff Dillon]

@Sam Chak Luckily I am a jet turbine expert, I've been working with these turbines for years. I can assure you that my description is accurate, I've measured the thrust response personally many times. The thrust is directly proportional to the turbine RPM. You can see videos on my website at https://www.purerockets.com. Totally understood about nonlinear systems now. If the behavior changes over time, it becomes a nonlinear system. So how would I change the Simulink plant model from the current (linear) transfer function to a nonlinear representation? Thanks for your time on this!

[Jeff Dillon]

This is a diagram of the turbine thrust response, given a step input command of 100% throttle

[Paul]

@Jeff Dillon, can you add at least one more curve to the plot showing the thrust response to a different step input of throttle command, like maybe 50%

@Sam Chak, linearity is a property of the input-output characteristic of the system, regardless of the representation of the system (ODE, etc.). Maybe you meant it's not feasible to test for linearity using a finite number of simulated input-output responses?

[Jeff Dillon]

@Paul If you give a 50% throttle command as a step input, it will take 2.5 seconds to reach 1/2 throttle. The graph I shared previously is accurate across all inputs. If you are already at 1/2 throttle for example, going to 100% will also take 2.5 seconds.

[Sam Chak]

Hi @Paul

Other than the input-output characteristic of the jet turbine, I have also asked if @Jeff Dillon can identify if the actuator has a nominal model like [linear] or [nonlinear]. I believe that the jet turbine itself has an internal control mechanism that regulates the throttle so that it can produce the output linearly with time t as shown.

At the moment, I don't know about model for the throttle dynamics, but I guess that the ODE is nonlinear. The following is not a model for the jet turbine, but an example of a nonlinear model (mathematically made, does not exist in real world) for a system that produces the 100% output at approximately 5 seconds for the input . The input-output characteristic is nonlinear in my example.

tspan  = linspace(0, 10, 10001);
x0     = [0 0.001];
[t, x] = ode45(@odefcn, tspan, x0);

% Plot
plot(t, x(:,1), 'linewidth', 1.5), 
grid on, xlabel('t, (sec)'), ylabel('Throttle')

% Dynamics of actuator
function xdot = odefcn(t, x)
    xdot    = zeros(2, 1);              
    u       = 1.0;                      % throttle percentage
    p1      = 2.8488;                   % tuning parameter
    p2      = 9.2/(p1*u);               
    xdot(1) = x(2);
    xdot(2) = p2*(1 - 2*x(1)/u)*x(2);   % nonlinear ODE
end

[Jeff Dillon]

I'm looking for an update to the Simulink plant model, not MATLAB code. I'm using the Simulink model in real-time with the actual hardware. There is no regulator inside the turbine. It accepts direct PWM input just like the motor on a quadcopter. A quadcopter motor also has a slight response delay, this is no different. A quadcopter motor has to spin the propeller and the delay is caused by the inertia of the propeller. It takes about about a 1/4 of a second for the quad rotor system to spin up to full throttle. The jet turbine is no different, you are making this unnecessarily difficult. When you provide a full throttle PWM input to the jet engine, it instructs the fuel pump to send more fuel. It then takes 5 seconds for the internal RPM of the jet turbine to reach full throttle. There is no regulator, it's just physics. For the turbine, the internal compressor spins at 10,000 RPM at zero throttle, and 140,000 RPM at full throttle when provided sufficient fuel from the fuel pump. If possible, please provide a corresponding update to the Simulink plant model, which is currently using a transfer function.

[Sam Chak]

Hi @Jeff Dillon

I apologize for any frustration you may be experiencing as a jet turbine expert regarding this modeling matter. As someone with a background in viewing time-dependent phenomena as dynamics, I tend to describe such phenomena using differential equations. These mathematical equations have been implemented in MATLAB and Simulink throughout my years of experience. If the use of differential equations is not the approach you wish to take in implementing the Simulink model to generate the desired throttle outputs, I suggest consulting with @Paul.

% Samples of throttle outputs

[Jeff Dillon]

@Sam Chak Sorry what do you mean "is not the approach I wish to take"? What approach are you suggesting? How do I use the differential equations that you provided in Simulink? Can you show the updated Simulink diagram that is using the differential equations in the plant model? I think I'm missing something here. How did you generate the images you just shared, are they outputs of the Simulink model?

[Sam Chak]

Hi @Jeff Dillon

What do you mean "is not the approach I wish to take"?

The ODE approach mentioned in my comment to @Paul is simply an example of a model for a hypothetical nonlinear actuator. It is not intended to be a model for the Jet Turbine. I apologize if you interpreted it as unnecessarily difficult for unknown reasons. That was not my intention.

What approach are you suggesting?

I don't want to make the ODE approach appear complicated for you. That's why suggest consulting with @Paul.

How do I use the differential equations that you provided in Simulink?

The ODEs were never intended for the Jet Turbine itself. They were used to demonstrate to @Paul my attempt at modeling the dynamics and generating the input-output characteristics. Therefore, they are not applicable to your Simulink modeling.

Can you show the updated Simulink diagram that is using the differential equations in the plant model?

I apologize for the confusion. I didn't save the Simulink file 9 months ago. However, if you are interested in working with nonlinear ODEs (without using transfer functions), I recommend checking out this video:

https://www.youtube.com/watch?v=A_AZd09Rml8

How did you generate the images you just shared?

They are plotted as a time-dependent piecewise algebraic function, without involving derivatives or ODEs. It is a simple math formula that describes the saturation characteristics.

Are they outputs of the Simulink model?

No, they are directly plotted in MATLAB. However, I assume these are the step response plots that @Paul requested from you. I hope they help Paul guide you in the next steps of the non-ODE modeling process in Simulink.

[Paul]

Hi Sam,

I don't really have any more to offer on this topic than you have already, which has all been directly on point, and certainly not with the information provided. I was just trying to make a minor point about the definition of linearity, which wasn't really relevant to the modeling discussion, and elicit more information about the system in question, which I thought might be relevant to the discussion.

[Jeff Dillon]

This problem doesn't seem much different than in a typical motor control system where the motor doesn't respond instantly to control input but instead follows a ramp

Using ChatGPT:

Me: I have a Simulink model where a jet turbine is mounted vertically and used to balance a VTOL system. The model is working using a PID controller but assumes that the jet turbine responds to control input instantly. In reality, the turbine takes 5 seconds to go from zero to full throttle, linearly. How would I modify the plant model which currently uses a transfer function to accommodate for this response delay?

ChatGPT: To modify your plant model in Simulink to accommodate for the response delay of the jet turbine, you can introduce a time delay block in the control loop. This time delay block will simulate the 5-second delay before the turbine responds to the control input. Here's how you can do it:

1. Open your Simulink model.
2. Locate the transfer function block that represents the jet turbine dynamics. This block should be connected to the PID controller.
3. Remove the existing transfer function block representing the turbine dynamics.
4. Add a time delay block in its place. You can find this block in the Simulink Library Browser under the "Discontinuities" category.
5. Double-click on the time delay block to open its parameters dialog box.
6. Set the time delay value to 5 seconds. This value represents the delay between the control input and the response of the turbine.
7. Connect the output of the time delay block to the rest of the model, as it was connected before with the transfer function block.
8. Save your Simulink model and run the simulation to observe the effect of the time delay on the system's response.

By adding the time delay block in the control loop, you are accounting for the actual response time of the jet turbine. This delay will introduce a more realistic behavior into your model, considering the 5-second time it takes for the turbine to go from zero to full throttle.

Me: Well, it's not really a 5 second delay before the turbine responds. It responds quickly, but takes 5 seconds to achieve full throttle. The thrust begins to increase immediately upon control input, but takes time to ramp up.

ChatGPT: I apologize for the misunderstanding. If the turbine responds quickly but takes 5 seconds to achieve full throttle, you can modify your plant model accordingly. Instead of using a simple time delay block, you can use a ramp-up function to simulate the gradual increase in thrust. Here's how you can do it:

1. Remove the time delay block that was added previously (if you have already added it).
2. Add a ramp-up function to your Simulink model. You can use a combination of blocks to achieve this effect. Here's an example of how you can construct the ramp-up function:a. Use a Unit Delay block (found in the "Discrete" category in the Simulink Library Browser) to hold the previous control input value.b. Subtract the previous control input from the current control input to obtain the change in input.c. Use an Integrator block (found in the "Continuous" category) to integrate the change in input over time.d. Multiply the output of the integrator by a gain value to control the rate of ramp-up.e. Add the output of the ramp-up block to the initial control input.
3. Adjust the gain value in the ramp-up block to control the rate at which the thrust ramps up. You may need to tune this value to match the desired time it takes for the turbine to achieve full throttle.
4. Connect the output of the ramp-up block to the rest of your Simulink model.

By using this ramp-up function, you are simulating the gradual increase in thrust over time, which aligns with the behavior of the turbine taking 5 seconds to achieve full throttle

[Jeff Dillon]

@Paul @Sam Chak Does this help? The directions are a bit over my head as I'm still learning Simulink, might either of you be able to follow up?

[Jeff Dillon]

How does the Simulink subsystem above fit into the plant model? Please see attached.

[Paul]

Out of curiosity, how does the turbine respsond to throttle inputs that are not step changes? For example, what is the response to a ramp input, or a small amplitude sine wave input, and how does that change with the frequency of the sine wave? Or are step changes the only type of throttle input relevant to this type of system?

[Jeff Dillon]

Yes, only step changes via PWM will be expected. Discreet PWM inputs. Examples are:

• From 0 to 10%
• 15% to 20%
• 75% to 74%
• 0% to 100%
• 100% to 0%
• ...

In all these cases, the first response graph that I provided would apply.

Were you able to open the Simulink diagram I provided? Any chance you could incorporate the subsystem that Sam provided? Will modifying the behavior of the plant model as we are discussing, also affect the controller model? Because I don't really "need" to model the plant, I'm ready to use actual hardware. I'm using Simulink in real time via ROS to control the hardware. That's all working. I just don't want this thing to fly off into space! In reality, it's confined to a test tower so there's no chance of damage or fly-away. There is no closed loop in this Simulink project I'm using in the video, only a StateFlow model with open loop control. But it's a start! https://www.youtube.com/watch?v=UM_sCuoT8jY

[Jeff Dillon]

@Sam Chak and @Paul Any updates? I'm actually changing the experiment a little. If you see in the video, I also have electric ducted fans for fast response during takeoff and landing. They provide 7 pounds of thrust each, there's 4 of them. What I plan to do is just use the turbine at a steady throttle just to offset the weight a bit, and then use the fans to control the altitude. Their response is fast enough to assume instantaneous behavior, so the previous Simulink model that you provided should work, with adjusted PID params. This test is closer to what the actual flight behavior will be anyway.

[Jeff Dillon]

@Sam Chak Your reply (including others, now deleted) are not relevant. I've explained in good detail regarding the input/response behavior of a jet engine. You don't need to know how it works for the same reason you don't need to know how a battery works at the electro-chemical level in order to model it's behavior. I'm not looking for a ramp controller. I simply stated that the throttle response of the jet turbine is slow. It's exactly like turning on a heater (with a step input) and the heater takes 7 seconds to warm up to full temperature, linearly. Or when you step on the gas in an automobile, it takes 7 seconds to achieve full speed. I need a PID controller that ACCOMODATES this response behavior, not induces it. And you ignored my last comment regarding the use of EDF's and apparently did not view the video I shared, the problem now allows instantaneous response using electric fans, similar to a quadcopter controller. The jet engine can be completely ignored, it's only used in steady state to offset the majority of the vehicle mass with the remainder of the control done by the electric fans. Google for "EDF electric ducted fan"

[Jeff Dillon]

@Sam Chak I get that you may be a little embarrassed on your misunderstanding of the problem statement and is why you removed your posts, hopefully it's clear now! Is there anyone else at Mathworks that you can ask? You had presented a non-linear MATLAB method previously that would likely apply, however you were not able to show it in the context of an updated Simulink model. It's ok. FYI, I'm following the course here http://apmonitor.com/pdc/index.php/Main/ArduinoTemperatureControl on Control Systems, hopefully I'll be able to answer my own questions in due time! Thanks for your help so far.