This example shows how to adjust prediction and control horizons at run-time to evaluate controller performance without recreating the controller object or regenerating the code.
Prediction and control horizons, together with controller sample time, are determined typically before other MPC settings such as constraints and weights are designed.
There are certain guidelines to help choose the sample time
Ts, prediction horizon
p, and control horizon
m. For example, assume you want to determine how far the controller should look into the future. In theory, prediction time should be long enough to capture the dominant dynamic behavior of the plant but not any longer so as to avoid wasting resources used in computation. In practice, you often start with a small value and gradually increase it to see how control performance improves. When it plateaus, stop.
Control horizon determines how many decision variables MPC uses in optimization. If the value is too small, you don't have enough degrees of freedom to achieve a satisfactory performance. On the other hand, if the value is too large, both computation load and memory footprint increase significantly with little performance improvement. Therefore, it is another place you want to try different values and compare the results.
In this example, we demonstrate how to adjust prediction and control horizons of an MPC Controller block using its inports and compare control performance after multiple runs of simulation without recreating MPC controller object used by the block. If the block is running on an embedded system, you can adjust the horizons in real-time too, without regenerating and redeploying the code.
To run this example, Simulink® and Simulink Control Design™ are required.
if ~mpcchecktoolboxinstalled('simulink') disp('Simulink is required to run this example.') return end if ~mpcchecktoolboxinstalled('slcontrol') disp('Simulink Control Design is required to run this example.') return end
The single-input-single-output nonlinear plant is implemented in Simulink model
mpc_nloffsets. At the nominal operating point, the plant is at steady state with output of -0.5.
plant_mdl = 'mpc_nloffsets';
operspec command from Simulink Control Design to create an operating point specification object with the desired output value fixed at steady state.
op = operspec(plant_mdl); op.Outputs.Known = true; op.Outputs.y = -0.5;
findop command from Simulink Control Design to obtain the nominal operating point.
[op_point, op_report] = findop(plant_mdl,op);
Operating point search report: --------------------------------- opreport = Operating point search report for the Model mpc_nloffsets. (Time-Varying Components Evaluated at time t=0) Operating point specifications were successfully met. States: ---------- <strong>Min</strong> <strong>x</strong> <strong>Max</strong> <strong>dxMin</strong> <strong>dx</strong> <strong>dxMax</strong> <strong>____</strong> <strong>_______</strong> <strong>___</strong> <strong>_____</strong> <strong>___________</strong> <strong>_____</strong> (1.) mpc_nloffsets/Integrator -Inf 0.59453 Inf 0 1.0258e-13 0 (2.) mpc_nloffsets/Integrator2 -Inf 2.1891 Inf 0 -1.0989e-09 0 Inputs: ---------- <strong>Min</strong> <strong>u</strong> <strong>Max</strong> <strong>____</strong> <strong>_______</strong> <strong>___</strong> (1.) mpc_nloffsets/In1 -Inf -1.1806 Inf Outputs: ---------- <strong>Min</strong> <strong>y</strong> <strong>Max</strong> <strong>____</strong> <strong>____</strong> <strong>____</strong> (1.) mpc_nloffsets/Out1 -0.5 -0.5 -0.5
linearize command from Simulink Control Design to linearize the plant at the nominal operating condition.
plant = linearize(plant_mdl, op_point);
Obtain nominal plant states, output and input.
x0 = [op_report.States(1).x;op_report.States(2).x]; y0 = op_report.Outputs.y; u0 = op_report.Inputs.u;
The linearized plant is underdamped second order system. Using the
damp command, we can find out the dominant time constant of the plant, which is about 1.7 seconds.
Pole Damping Frequency Time Constant (rad/seconds) (seconds) -5.95e-01 + 1.84e+00i 3.07e-01 1.94e+00 1.68e+00 -5.95e-01 - 1.84e+00i 3.07e-01 1.94e+00 1.68e+00
A simple guideline recommends that the prediction time should at least cover the dominant time constant (1.7 seconds) and control horizon is 10%~20% of the prediction horizon. Therefore, if we choose sample time of 0.1, the prediction horizon should be around 17. This gives us a starting point to choose the default horizons
Ts = 0.1; p = 20; m = 4; mpcobj = mpc(plant,Ts,p,m);
-->The "Weights.ManipulatedVariables" property of "mpc" object is empty. Assuming default 0.00000. -->The "Weights.ManipulatedVariablesRate" property of "mpc" object is empty. Assuming default 0.10000. -->The "Weights.OutputVariables" property of "mpc" object is empty. Assuming default 1.00000.
Set nominal values in the controller.
mpcobj.Model.Nominal = struct('X', x0, 'U', u0, 'Y', y0);
Set MV constraint.
mpcobj.MV.Max = 2; mpcobj.MV.Min = -2;
Since there is little noise in the plant, we reduce the noise model gain to make the default Kalman filter more aggressive.
mpcobj.Model.Noise = 0.1;
mpc_onlineHorizons model implements the closed-loop control system. Our goal to track a -0.2 step change in the reference signal with minimum overshoot. We also want the settling time to be less than 5 seconds.
r0 = -0.7; mdl = 'mpc_onlineHorizons'; open_system(mdl)
In the model, the MPC block has two inports where we can connect prediction horizon signal and control horizon signal. In each simulation, we vary the prediction horizon value (from 5 to 50) while keeping the control horizon at 4. We measure both the overshoot (%) and settling time (sec) from the saved simulation results. Note that the MPC controller object is not changed. Instead, the new horizon values are supplied as input signals at run-time.
p_choices = 5:5:50; set_param([mdl '/Control Horizon'],'Value','4') for p = p_choices set_param([mdl '/Prediction Horizon'],'Value',num2str(p)) sim(mdl,20) settling_timeP(p/5) = ... find((abs(y.signals.values-r0)<0.01)&(abs([0;diff(y.signals.values)])<0.001),1,'first')*Ts; if r0>y0 overshootP(p/5) = abs((max(y.signals.values)-r0)/r0)*100; else overshootP(p/5) = abs((min(y.signals.values)-r0)/r0)*100; end end figure subplot(2,1,1) plot(p_choices,overshootP,'*') xlabel('prediction horizon') ylabel('overshoot (%)') title('control horizon = 4') subplot(2,1,2) plot(p_choices,settling_timeP,'*') ylabel('settling time (sec)') xlabel('prediction horizon')
-->Converting model to discrete time. -->Assuming output disturbance added to measured output channel #1 is integrated white noise.
As the two plots show above, when prediction horizon increases from 5 to 15, the overshoot drops from 6% to 3% and settling time increases from 3 seconds to 4 seconds. After that, however, both overshoot and settling time remain more or less the same. In addition all the settling time values satisfy the upper bound of 5 seconds. Therefore, we choose the prediction horizon of 15, because it is the smallest value to achieve satisfactory performance by forming the smallest optimization problem.
After we choose the prediction horizon, we use the same setup to evaluate different control horizon choices. In each simulation, we vary the control horizon (from 1 to 10) while keeping the prediction horizon at 15.
c_choices = 1:10; set_param([mdl '/Prediction Horizon'],'Value','15') for c = c_choices set_param([mdl '/Control Horizon'],'Value',num2str(c)) sim(mdl,20) settling_timeC(c) = ... find((abs(y.signals.values-r0)<0.01)&(abs([0;diff(y.signals.values)])<0.001),1,'first')*Ts; if r0>y0 overshootC(c) = abs((max(y.signals.values)-r0)/r0)*100; else overshootC(c) = abs((min(y.signals.values)-r0)/r0)*100; end end figure subplot(2,1,1) plot(c_choices,overshootC,'*') xlabel('control horizon') ylabel('overshoot (%)') title('prediction horizon = 15') subplot(2,1,2) plot(c_choices,settling_timeC,'*') xlabel('control horizon') ylabel('settling time (sec)')
As the two plots show above, when control horizon increases from 1 to 3, the overshoot drops from 10% to 2%. After that, it increases back to 5% as control horizon grows from 4 to 10. The explanation is that when control horizon is 1, the controller doesn't have enough degrees of freedom to achieve reasonable response. When control horizon is 4 or beyond, the controller has more decision variables such that the first optimal move often becomes more aggressive and thus results in larger overshoot but shorter settling time. In this example, since the main control goal is to achieve minimum overshoot, we choose 3 as control horizon.
The model is simulated with prediction horizon = 15 and control horizon = 3. Recall that our original design choice is prediction horizon = 20 and control horizon = 4 based on a simple guideline, which is close to our final choice.
set_param([mdl '/Prediction Horizon'],'Value','15') set_param([mdl '/Control Horizon'],'Value','3') open_system([mdl '/Input']) open_system([mdl '/Output']) sim(mdl)
The major benefit of using run-time prediction and control horizon inports in MPC and Adaptive MPC blocks is that you can evaluate and adjust controller performance in real-time without regenerating code and re-deploying it to the target system. This feature is very helpful at the prototyping stage.
To use run-time horizon adjustment in real time, the target system must support dynamic memory allocation because as horizons change, the sizes of all that matrices used to construct the optimization problem change at run-time as well.
You also need to specify the maximum prediction horizon in the block dialog to define the upper bound of the sizes of these matrices. Therefore, the memory footprint would be large. After finding the best horizon choices, it is recommended to disable the feature to have efficient code generation with fixed-size data.