About Linear Parameter Varying (LPV) Models

In many situations the nonlinear dynamics of a system need to be approximated using simpler linear systems. A single linear system provides a reasonable model for behavior limited to a small neighborhood around an operating point of the nonlinear system. When the nonlinear behavior needs to be approximated over a range of operating conditions, we can use an array of linear models that are interconnected by suitable interpolation rules. Such a model is called an LPV model.

To generate an LPV model, the nonlinear model is trimmed and linearized over a grid of operating points. For this purpose, the operating space is parameterized by a small number of scheduling parameters. These parameters are often a subset of the inputs, states, and output variables of the nonlinear system. An important consideration in the creation of LPV models is the identification of a scheduling parameter set and selection of a range of parameter values at which to linearize the model.

We illustrate this approach for approximating the pitch dynamics of an airframe.

Pitch Dynamics of an Airframe

Consider a three-degree-of-freedom model of the pitch axis dynamics of an airframe. The states are the Earth coordinates , the body coordinates , the pitch angle , and the pitch rate . Figure 1 summarizes the relationship between the inertial and body frames, the flight path angle , the incidence angle , and the pitch angle .

Figure 1: Airframe dynamics.

The airframe dynamics are nonlinear and the aerodynamic forces and moments depend on speed and incidence . The model scdairframeTRIM describes these dynamics.

open_system('scdairframeTRIM')

Batch Linearization Across the Flight Envelope

Use the speed and the incidence angle as scheduling parameters; that is, trim the airframe model over a grid of and values. Note that these are two of the five outputs of the scdairframeTRIM model.

Assume that the incidence varies between -20 and 20 degrees and that the speed varies between 700 and 1400 m/s. Use a 15-by-12 grid of linearly spaced pairs for scheduling:

nA = 15;   % number of alpha values
nV = 12;   % number of V values
alphaRange = linspace(-20,20,nA)*pi/180;
VRange = linspace(700,1400,nV);
[alpha,V] = ndgrid(alphaRange, VRange);

For each flight condition , linearize the airframe dynamics at trim (zero normal acceleration and pitching moment). This requires computing the elevator deflection and pitch rate that result in steady and .

Use operspec to specify the trim condition, use findop to compute the trim values of and , and linearize the airframe dynamics for the resulting operating point. See the "Trimming and Linearizing an Airframe" example for details.

The body coordinates, , are known states for trimming. Therefore, you need to provide appropriate values for them, which you can specify explicitly. However, in this example, let the model derive these known values based on each pair. For each flight condition , update the values in the model and create an operating point specification. Repeat these steps for all 180 flight conditions.

clear op report
for ct = 1:nA*nV
   alpha_ini = alpha(ct);      % Incidence [rad]
   v_ini = V(ct);              % Speed [m/s]

   % Specify trim condition
   opspec(ct) = operspec('scdairframeTRIM');

   % Xe,Ze: known, not steady.
   opspec(ct).States(1).Known = [1;1];
   opspec(ct).States(1).SteadyState = [0;0];

   % u,w: known, w steady
   opspec(ct).States(3).Known = [1 1];
   opspec(ct).States(3).SteadyState = [0 1];

   % theta: known, not steady
   opspec(ct).States(2).Known = 1;
   opspec(ct).States(2).SteadyState = 0;

   % q: unknown, steady
   opspec(ct).States(4).Known = 0;
   opspec(ct).States(4).SteadyState = 1;

end
opspec = reshape(opspec, [nA nV]);

Trim the model for all of the specified operating point specifications.

Options = findopOptions('DisplayReport','off', ...
	'OptimizerType','lsqnonlin');
Options.OptimizationOptions.Algorithm = 'trust-region-reflective';
[op, report] = findop('scdairframeTRIM',opspec,Options);

The op array contains the operating points found by findop that will be used for linearization. The report array contains a record of input, output, and state values at each point.

Specify linearization inputs and outputs.

io = [linio('scdairframeTRIM/delta',1,'in');...        % delta
   linio('scdairframeTRIM/Airframe Model',1,'out');... % alpha
   linio('scdairframeTRIM/Airframe Model',2,'out');... % V
   linio('scdairframeTRIM/Airframe Model',3,'out');... % q
   linio('scdairframeTRIM/Airframe Model',4,'out');... % az
   linio('scdairframeTRIM/Airframe Model',5,'out')];   % gamma

Batch-linearize the model at the trim conditions. Store linearization offset information in the info structure.

[G,~,info] = linearize('scdairframeTRIM',op,io, ...
	linearizeOptions('StoreOffsets',true));
G = reshape(G,[nA nV]);
G.u = 'delta';
G.y = {'alpha','V','q','az','gamma'};
G.SamplingGrid = struct('alpha',alpha,'V',V);

G is a 15-by-12 array of linearized plant models at the 180 flight conditions . The plant dynamics vary substantially across the flight envelope, including scheduling locations where the local dynamics are unstable.

bodemag(G(3:5,:,:,:))
title('Variations in airframe dynamics')

The LPV System Block

The LPV System block in the Control System Toolbox™ block library facilitates simulation of linear parameter varying systems. The primary data required by the block is the state-space system array G that was generated by batch linearization. We augment this with information about the input/output, state, and state derivative offsets from the info structure.

Extract the offset information.

offsets = getOffsetsForLPV(info);
xOffset = offsets.x;
yOffset = offsets.y;
uOffset = offsets.u;
dxOffset = offsets.dx;

LPV Model Simulation

Open the system scdairframeLPV, which contains an LPV System block that has been configured based on linear system array G and the various offsets.

open_system('scdairframeLPV')

An input signal was prepared based on a desired trajectory of the airframe. This signal u and corresponding time vector t are saved in the scdairframeLPVsimdata.mat file. Specify the initial conditions for simulation.

alpha_ini = 0;
v_ini = 700;
x0 = [0; 700; 0; 0];
sim('scdairframeLPV')

The simulation shows good emulation of the airframe response by the LPV system. We chose a very fine gridding of scheduling space leading to a large number (180) of linear models. Large array sizes can increase implementation costs. However, the advantage of LPV representations is that we can adjust the scheduling grid (and hence the number of linear systems in the array) based on:

The former information helps reduce the range for the scheduling variables. The latter helps pick an optimal resolution (spacing) of samples in the scheduling space.

Let us plot the actual trajectory of scheduling variables in the previous simulation against the backdrop of gridded scheduling space. The outputs were logged via their scopes (contained inside the Compare Responses block of scdairframeLPV).

Stable = false(nA,nV);
for ct = 1:nA*nV
   Stable(ct ) = isstable(G(:,:,ct));
end
alpha_trajectory = Alpha_V_Data.signals(1).values(:,1);
V_trajectory = Alpha_V_Data.signals(2).values(:,1);

plot(alpha(Stable)*180/pi,V(Stable),'g.',...
   alpha(~Stable)*180/pi,V(~Stable),'k.',...
   alpha_trajectory,V_trajectory,'r.')
title('Trajectory of scheduling variables')
xlabel('\alpha'); ylabel('V')
legend('Stable locations','Unstable locations','Actual trajectory')

The trajectory traced during simulation is shown in red. Note that it traverses both the stable and unstable regions of the scheduling space. Suppose you want to implement this model on a target hardware for input profiles similar to the one used for simulation above, while using the least amount of memory. The simulation suggests that the trajectory mainly stays in the 890 to 1200 m/s range of velocities and -15 to 12 degree range of incidence angle. Furthermore, you can explore increasing the spacing between the sampling points. Suppose you use only every third sample along the dimension and every second sample along the dimension. The reduced system array meeting these constraints can be extracted from G as follows:

I1 = find(alphaRange>=-15*pi/180 & alphaRange<=12*pi/180);
I2 = find(VRange>=890 & VRange<=1200);
I1 = I1(1:2:end);
I2 = I2(1:3:end);

Gr = G(:,:,I1,I2);
size(Gr)

5x2 array of state-space models.
Each model has 5 outputs, 1 inputs, and 4 states.

The new sampling grid, Gr, has a more economical size of 5-by-2. Simulate the reduced model and check its fidelity in reproducing the original behavior.

Change directory to a writable directory since model would need to be recompiled

cwd = pwd;
cd(tempdir)
lpvblk = 'scdairframeLPV/LPV System';
set_param(lpvblk,...
   'sys','Gr',...
   'uOffset','uOffset(:,:,I1,I2)',...
   'yOffset','yOffset(:,:,I1,I2)',...
   'xOffset','xOffset(:,:,I1,I2)',...
   'dxOffset','dxOffset(:,:,I1,I2)')
sim('scdairframeLPV')
cd(cwd)

There is no significant reduction in overlap between the response of the original model and its LPV proxy.

The LPV model can serve as a proxy for the original system in situations where faster simulations are required. The linear systems used by the LPV model may also be obtained by system identification techniques (with additional care required to maintain state consistency across the array). The LPV model can provide a good surrogate for initializing simulink design optimization problems and performing fast hardware-in-loop simulations.