LFRT version 1.3 for Matlab: introduction

This tutorial is not the best way to discover the LFRT toolbox because the modelling example that is treated concentrates on potential problems. Simpler illustrative examples can be found in the manual. They are also available in the sub-directory lfrtemo of the installation directory.

In this page, we consider a classical missile model, we shall treat the following items

• Modelling a simplified missile
• Comments on modelling
• Scheduled feedback design (LFT form) and analysis
• Comments on feedback gains in LFT form
• Alternative scheduled feedback design

Example 1 : Modelling a simplified missile

• Non-linear equations of the missile
• LFT representation of the continuum of linearized models
• Models for feedback synthesis and for analysis
• LFT representation of these models

1- Non-linear model: Using standard notation (α angle of incidence, q pitch rate, Ma Mach number, δp tail plane deflection), the nonlinear longitudinal equation is

The input is δp. There are two measurements, the load factor (η) and the pitch rate (q)

The numerical values are defined here. Note that only α (denoted Al), q, Ma, δp (denoted dp) do not have a numerical value (other numerical values are initialized by invoking the LFRT script  missiledata, this file is in the directory lfrdemo).

2- LFT representation of the continuum of linearized models: In the first part of this script we use the symbolic toolbox in order to compute linearized models. Note that the LFRT toolbox can also be used in a similar way (for differentiation using difflfr, for symbolic evaluation (eval) using pluglfr). But the use of the symbolic toolbox permits us to take advantage of the tree decomposition for realization (symtreed). Using symtreed and then minlfr as proposed here, is often the most efficient approach to low order modelling.

% Load numerical data
   missiledata

% Define symbolic objects
   syms Al q Ma dp

% Build differential equations
   Cz = z3*Al^3 + z2*Al^2 + z1*( 2 -(1/3)*Ma)*Al + z0*dp;
   Cm = m3*Al^3 + m2*Al^2 + m1*(-7 +(8/3)*Ma)*Al + m0*dp;

   A1 = q+K1*Ma*Cz*(1-Al^2/2);%+Al^4/24);
   A2 = K2*Ma^2*Cm;
   C1 = K3*Ma^2*Cz;

% Differentiate for obtaining linearized models
   A = [diff(A1,'Al') diff(A1,'q');diff(A2,'Al') diff(A2,'q')];
   B = [diff(A1,'dp');diff(A2,'dp')];

   C = [0 1;diff(C1,'Al') diff(C1,'q') ];
   D = [0 ; diff(C1,'dp') ];

% Symbolic form of A,B,C,D
   dp = -(m3*Al^3 + m2*Al^2 + m1*(-7 +(8/3)*Ma)*Al)/m0;
   A = eval(A); B = eval(B); C = eval(C); D = eval(D);

An lfr-realization of the system matrix [A B;C D] is performed using the function symtreed (tree decomposition), and then, the input/output corresponding form is computed using abcd2lfr. The result is further reduced using minlfr.

% Transformation to lfr-object
   ABCD = symtreed([A B;C D],[Al Ma],[0.17 3]);
   sys = abcd2lfr(ABCD,2);
   sys = minlfr(sys,1000*eps);

Now, the symbolic object representing the continuum of linearized models is an lfr-object (sys).

Normalisation: The angle of incidence is assumed to vary in [0  0.24], and the Mach number in [2  4]. The mean values are respectively 0.17 and 3. When we invoked symtreed, the argument [0.17 3] defined a translation in the parameter space so that the origin of variation is now the point [0.17  3]. It remains to normalize the variations between -1 and +1 by invoking normlfr as follows.

% Normalisation 
   sys = normlfr(sys,[-0.17  -1],[0.17  1]);

% Actuators added at the system input
   sys = sys*ss(tf([omegaa^2],[1 2*kia*omegaa omegaa^2]));

% Al and Ma defined as lfr-objects
   lfrs Al Ma

3- Models for feedback synthesis and for analysis: The following figure shows the controller architecture.

4- LFR representation of the models for synthesis and analysis: This part of the script illustrates the concatenation of lfr-objects.

% General model (this equation)
   mat1 = [1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 lfrS];
   mat2 = [ 1 0; [0;0;0] [1 0;0 1;0 -1]*sys];
   sys_gene  = mat1*mat2;

% Model for synthesis (from second input to measurements 2,3,and 4)
   sys_syn = diag([0 1 1 1])*sys_gene*[0;1];

% Model for analysis (assume that the feedback is klfr)
   % close the feedback loop
   sys_ana = feedback(sys_gene,klfr,1,1,3);
   % consider from first input to outputs 1 and 3
   sys_ana = diag([1 0 1 0])*sys_ana*[1;0];

   ssS = ss(0,1,1,0);
   distlfr(lfrS,ssS)

Comments on modelling

• When nothing works, or results are far from expected (e.g., the tree decomposition leads to a huge lfr-object), in most cases it is a problem of scaling. For example, for aircraft modelling, small values (radians) are mixed to huge values (mass in kg). Before using any function of the toolbox, replace natural units by artificial units so that parameter variations remain roughly around 1 or 10.
• There is some ambiguity behind the word "parameter". Let us define a parameter x as follows:

    lfrs x  

  In fact x, is not a parameter but an lfr-object of the form x = 1 δx (1 - 0 δx)^-1 1 + 0. Let us normalize x so that when δx varies in [-1  1], x varies in [1  4]:

     x = normlfr(x,1,4)  

  now x is an lfr-object of the form x = 1 δx (1 - 0 δx)^-1 1.5 + 2.5.

  The actual parameter is δx that coincide with x in the first case. In the second case, x and δx are different. However, we do not distinguish between both concepts, both are called "parameters" assuming that which one is considered is always clear from the context.

• Normalization of parameters is a source of misunderstanding. A normalization consists of two steps (treated often together) that are, first, a translation of the origin in the parameter space, second, a normalization between -1 and +1 (see above example). We need the first step in order to avoid division by zero. For example, consider a parameter x that belongs to [0   4]. The following (natural) script will not work

    lfrs x
    y = 3 + 1/x;
    y = normlfr(y,0,4);

  The problem using this script is that the computation of 1/x induces a division by zero. However, it is possible to perform the above computation by setting the mean value to 2 before computing 1/x (translation of the origin from 0 to 2). So, the correct version of the above script is

    lfrs x 2
    y = 3 + 1/x;
    y = normlfr(y,-2,2);

  The same problem was encountered in the missile example when symtreed was invoked (see here).

  This problem of division by zero will be treated automatically in version 2 of this toolbox (release end of 2004)

• There are some limitations due to the fact that the parameters don't have a name as in the Symbolic Toolbox but are recognized by ordering them. For example:

    lfrs a b c
    lfrs x y z
    M1 = a*(y^2+z^2);
    M2 = x*(b^2+c^2);
    distlfrs(M1,M2)
      ans =
          0.

  shows that M1 and M2 are the same object, similarly, a is x, b is y and c is z.

  This problem was encountered here: After invoking symtreed the parameters exist inside the lfr-object created by this function but are not available outside. That is why we re-define them in the workspace using lfrs.

  This problem will disappear with version 2 of this toolbox (release end of 2004) because all parameters will be recognized by their name rather than by ordering.

Example 2 : Scheduled feedback design (LFT form) and analysis

Let us consider the system sys_syn defined earlier. This system has three dominant modes, as there are three measurements we can use directly fb_sched in order to assign three poles by a feedback gain in LFT form.

The technique behind fb_sched is eigenstructure assignment (see here). Assuming that a feedback gain is a function K of (A,B,C,D) denoted K(A,B,C,D), the LFT gain computed by this function is K(Δ) = K(A(Δ),B(Δ),C(Δ),D(Δ)) in which A(Δ), B(Δ) and so on stand for the LFT forms of A and B. The computation of K(A(Δ),B(Δ),C(Δ),D(Δ)) is feasible because the function K(A,B,C,D) only involves matrix operations that have an LFT counterpart. (Note that alternative design techniques, e.g. techniques using eigenvalue computation (like LQG) cannot work because eigenvalues or LFTs are not LFTs.)

• A pole will be assigned to -10 for all systems and
• a pair of poles will be considered as an LFT: -19+4*Ma+(19+4*Al)*j (and conjugate) where Al and Ma (already normalized) vary from -1 to +1, it means that the pole with positive imaginary part will vary in a rectangle [-15  -23] + [15  23]*j (obviously this assignment is just for fun, it has no physical meaning).

% Name given to the (normalized) parameters of the missile model
   lfrs Al Ma

% Scheduled feedback design
   pol = [-10 -19+4*Ma+(19+4*Al)*j];
   klfr = fb_sched(sys_syn,0,pol);

% Verification: closed-loop poles at some points of the flight domain
   [eig_fb(sys_syn,klfr,[0 0]) ...
   eig_fb(sys_syn,klfr,[1 1]) ...
   eig_fb(sys_syn,klfr,[-1 -1])]

      ans =

      1.0e+02 *

      -0.8156 + 0.4250i  -0.8575 + 0.7349i  -1.3477
      -0.8156 - 0.4250i  -0.8575 - 0.7349i  -0.2300 + 0.1500i
      -0.1900 + 0.1900i  -0.1500 + 0.2300i  -0.2300 - 0.1500i
      -0.1900 - 0.1900i  -0.1500 - 0.2300i  -0.1976
      -0.1000            -0.1000            -0.1000          

The feasibility of the on-line computation of such a gain is not guaranteed because it has the following form:

    K(Δ) = K12Δ (I - K11Δ)^-1K12 + K22

% Well-posedness analysis of the gain
   [rad_min,rad_max,pert] = wp_rad(klfr)

      rad_min =

      1.7129

      rad_max =

      Inf

The LFRT toolbox proposes several tools based on parameter gridding in order to analyse directly the properties of a feedback loop, when the system or the controller or both are lfr-objects. First is presented a tool for pole maps.

% Parameter gridding
   grid = makegrid({1,-1,1,5},{2,-1,1,5});

% Symbols
   symb = symbplot(grid,{1,'dis','Al'},{2,'col','Ma'});

% Pole map
   figure
   ppole(sys_syn,klfr,grid,symb);
   axis([-80 5 -5 80]);

The LFRT toolbox proposes in particular an interface with the standard Matlab function ltiview. Here we illustrate step responses, but it is also possible to draw Bode, Nyquist and Nichols loci with this tool.

% Model for step response analysis (see here)
   sys_ana = feedback(sys_gene,klfr,1,1,3);
   sys_ana = diag([1 0 1 0])*sys_ana*[1;0];

% Symbols
   symb2 = symbplot(grid,{1,'con','Al'},{2,'col','Ma'});

   % Step responses
   lfrview('step',sys_ana,grid,symb2);

Comments on feedback gains in LFT form

• When there are less measurements than dominant modes.

  In this case, we suggest to add an observer of order equal to the difference between the number of dominant modes and the number of measurements. Note that all dominant modes must be treated because they are well controllable, so, sensitive to the treatment of other dominant modes.

  Observers are not used for estimating the state but for providing additional degrees of freedom to the designer. For this reason, the classical concept of minimum order observer is not relevant (using a first-order observer might be sufficient even if the so-called "minimum order" is much larger).

  The functions that permit us to proceed as above are ob_sched (observer design) addobs (connects the observer to the system) and obs2klfr (transforms to a non-observer-based equivalent scheduled controller). The non-LFT observer theory presented in the Robust Modal Control Toolbox (see function ob_gene) is sufficient for understanding the use of these three functions.

  The missile is treated with an observer (assuming tha q is not measured) in the file ex_8_2.m available in the sub-directory lfrtemo of the installation directory.

• Multi-inputs case.

  The missile example is single-input, therefore, for fixing all the degrees of freedom, it suffices to choose the closed-loop eigenvalues. In the multi-inputs case, choosing poles is not sufficient, so, some options for selecting the corresponding eigenvectors are proposed (see fb_sched, and for a comparison with the non-LFT case see: fb_prop).

• Order reduction of LFT gains.

  The order of LFT gains (size of the Δ matrix) might be very large even after using the order reduction function minlfr. This problem is often easily solved by considering an LFT representation specific to the feedback implementation. Example: to be written.

• Robustness and gain scheduling.

  In a system there are two kinds of varying parameters:

  • measured varying parameters <-> scheduling
  • uncertain parameters <-> robustness

  of course, it is necessary to treat simultaneously scheduling and robustness.

  The Robust Modal Control Toolbox was developed for treating robustness against real parameter variations. The key idea was to treat together a set of worst case models (detected sequentially) by a multimodel control design technique. The same idea can be applied for LFT gain design. In order to have enough degrees of freedom for assigning poles relative to a set of LFT models, we propose to use observers (see above).

  The key function for scheduling, but also for multimodel control design is fb_sched. In Example 2 it is used for scheduling, in Example 3 it is used for multimodel design. Of course, it can be used for both, simultameously.

• Discrete time implementation.

  To be written

Example 3 : Alternative scheduled feedback design

The LFT approach presented earlier is very powerful as it computes the feedback gain everywhere in the flight domain in one shot. However, LFT-gains obtained in that way are usually of large order (size of the Δ-matrix). In addition, there are on-line inversions (unless the discrete time form is implemented, in which case this problem becomes a problem of convergence of a fixed-point algorithm).

The alternative technique (Magni, 1999, Aerospace Science and Technology, vol 3(3)) presented here is less powerful as it consists of several steps before leading to a satisfactory scheduled law but the resulting gains have low order and do not require on-line inversion.

• A scheduling form of the feedback gain is chosen a priori, for example

     

• The system is augmented as follows (using the function lfr2bsys)

     

  because, it is clear that

     

• The original gain scheduling problem is now transformed to a problem of robustness because the feedback gain

     

  must be designed so that it is satisfactory for all augmented systems (we shall use the Robust Modal Control Toolbox for that purpose).
   
• Go back to original system by computing K(Δ) form (using the function bfb2lfr)

     

Choice of the interpolation formula and augmented system. We look for a feedback gain as follows K(Δ) = K0 + Al K1 + Al^2 K2 + Ma K3 + Ma^2 K4. Then, the corresponding augmented system is computed by invoking lfr2bsys. Let us consider the system sys_syn defined earlier.

% Name given to the (normalized) parameters of the missile model
     >> lfrs Al Ma

% Desired scheduling formula
     >> lfrex = [1 Al Al^2 Ma Ma^2];

% Corresponding augmented system
     >> sys_aug = lfr2bsys(sys_syn,lfrex)

Initial nominal feedback design. First is computed a simple gain (with K1 = K2 = K3 = K4 = 0) relative to the nominal system.

% STEP 1: Nominal gain
     >> k_aug_0 = fb_sched(sys_aug,0,[-10, -25+20j],{[0,0],[0,0]});

Transformation to LFT gain and analysis. Before using analysis tools it is necessary to transform this artificial gain into an LFT gain relative to the original system. For that purpose, the function bfb2lfr is invoked.

% Corresponding LFT gain
     >> klfr0 =  bfb2lfr(k_aug_0,lfrex);

% Pole map analysis
     >> grid = makegrid({1,-1,1,5},{2,-1,1,5});
     >> symb = symbplot(grid,{1,'dis','Al'},{2,'col','Ma'});
     >> figure
     >> ppole(sys_syn,klfr0,grid,symb);

Multimodel feedback design. We give below the result of our first trial. Better results with the relevant explanations can be found in (Döll et al, 2001, Control Engineering and Practice, vol 9, pp 1067-1078).

% pol is the set of assigned poles
% par defines the values of [Al Ma] for each assigned pole
     >> pol = [  -10,-25+20j,  -10,-30+30j,   -10,-30+30j,    -10,-10+10j];
     >> par = {[0,0],  [0,0],[1 1],  [1 1],[-1 1], [-1 1],[-1 -1],[-1 -1]};

     >> k_aug = fb_sched(sys_aug,0,pol,par);

fb_sched can also be used for mixing scheduling and multimodel control.

Transformation to LFT gain and analysis.

% Corresponding LFT gain
     >> klfr =  bfb2lfr(k_aug,lfrex);

% Pole map analysis
     >> figure
     >> ppole(sys_syn,klfr,grid,symb);

% Close the feedback loop
     >> sys_ana = feedback(sys_gene,klfr,1,1,3);

% Model for simulation
     >> sys_ana = [1 0 0 0;0 0 1 0]*sys_ana*[1;0];

% Step responses
     >> symb2 = symbplot(grid,{1,'con','Al'},{2,'col','Ma'});
     >> lfrview('step',sys_ana,grid,symb2);

µ-analysis. The system and the feedback are in LFT form. Therefore, it is straightforward to use µ-analysis to check the stability of the continuum of models:

% Upper and lower bounds
     >> frequ = logspace(0,3,100);
     >> [M,blk] = lfr2mua(sys_ana,frequ);
     >> bnds = mu(M,blk,'wlu');

% Plotting upper and lower bounds of mu
     >> figure
     >> semilogx(frequ,bnds(1:100,1),'r-'); hold on
     >> semilogx(frequ,bnds(1:100,2),'g.');

% Alternative lower bound (only if two dominant uncertainties)
     >> [Al_Ma,lmu,frequ] = lowermu(sys_ana)

Al_Ma =

   -0.7276    1.3096


lmu =

    0.7636


frequ =

     0

% Closed-loop eigenvalues of worst case
      >> eig(sys_ana,Al_Ma)

ans =

 -89.4742 +39.0204i
 -89.4742 -39.0204i
   0.0866
 -15.9457 +18.6739i
 -15.9457 -18.6739i

Improvement of this control law. The above technique is interesting for obtaining quickly a scheduled gain that is not too bad. For more ambitious design objectives it is better to use an automatic tuning technique (see function fb_tun) rather than a manual procedure based on trials and errors.