RMCT for Matlab: a tutorial

Before using the toolbox, it is recommended to read at least pages 109 to 115 of the manual. The RMCT Toolbox proposes three main tools for control design. These tools are

• fb_prop for proportional (constant gains) feedback design. Dynamic feedback design is also possible using this function if an observer has been connected to the system. Two functions combine observers and the use of fb_prop, these functions are dfb_proj and dfb_ins but in this tutorial, we shall only present the most basic tools.
• fb_tun for tuning a proportional feedback. Can be extended to observer-based dynamic feedback tuning. For using this function, the user must write a short script containing a loop.
• fb_dyn for computing a dynamic feedback that performs one or all the following objectives: multi-model eigenstructure assignment, multi-model "phase control", feedback gain structuring, frequency domain constraints.

• Example 1: elementary use of fb_prop - comments.
• Example 2: use of fb_prop with an observer - comments.
• Example 3: elementary use of fb_tun - comments.
• Example 4: use of fb_dyn for multimodel "phase control".
• Example 5: use of fb_dyn for multimodel eigenstructure assignment.
• Example 6: use of fb_dyn for controller order reduction - comments.

Example 1 : elementary use of fb_prop

The input arguments of this function are

• sys : linear continuous time system.
• k0 : a reference feedback. The distance of the feedback being computed to the reference feedback is minimized if assignment constraints do not fix all the degrees of freedom.
• pol_cl : closed-loop poles to be assigned, e.g. pol_cl = [-2,-4+4j,-5].
• key : option for eigenvector assignment ('m','p','v','z','n' or 'e'). If the same option is used for all assignments use key = 'p' otherwise key = 'pmpv' for 4 eigenvector assignments with respective options 'p','m','p','v'. See here for details.
• def_pb : matrix, each column of which defines the eigenvector to be assigned. The syntax for each column of this matrix depends on the corresponding option in key. See here for details.

The considered example consists of shifting with "minimum energy" 4 poles of a system. The considered system has 4 dominant poles and 4 measurements, so, the design procedure is very simple, it reduces to:

   sys = rcamdata('lat',0,1);
   pol_ol = [-0.23+0.59j -0.18 -1.30];
   pol_cl = [-0.60+0.60j -0.60 -1.30];
   K1 = fb_prop(sys,0,pol_cl,'m',pol_ol)

% Load a system with 6 states, 4 measurements and 2 inputs

   sys = rcamdata('lat',0,1);

% Compute open-loop eigenvalues

   eig(sys)

   ans =

     -1.3017
     -0.1837
     -0.2360 + 0.5954i
     -0.2360 - 0.5954i
     -6.6667
     -3.3333

% This system has 4 dominant poles, so, 4 poles must be treated
%    -0.23+0.59j  replaced by  -0.60+0.60j
%    -0.18        replaced by  -0.60
%    -1.30        preserved

   pol_ol = [-0.23+0.59j -0.18 -1.30]; % This is def_pb
   pol_cl = [-0.60+0.60j -0.60 -1.30];

% Assignment by invoking fb_prop. Pole shifting is performed with
% minimum energy option ('m')

   K1 = fb_prop(sys,0,pol_cl,'m',pol_ol)

% Results (eig_fb computes closed-loop eigenvalues)

   K1 =

       -0.7506    0.2549    0.8871    0.2580
       -0.7943    0.1637    2.0341    0.0377

   eig_fb(sys,K1)

   ans =

     -6.4311
     -2.4262
     -0.6000 + 0.6000i
     -0.6000 - 0.6000i
     -0.6000
     -1.3000

Additional comments on proportional feedback design.

• Alternative options: Alternative eigenvector selection options can be used, click here for details.
• Reference feedback case: The second input argument can be useful when we do not want to treat all dominant poles and have a reference feedback. For example if K1 (previously computed) is the reference feedback and if the eigenvalue assigned at -0.6+0.6j must be shifted to -0.7+0.7j we shall invoke fb_prop as follows:

  K = fb_prop(sys,K1,-0.7+0.7j,'m',-0.23+0.59j)

• Eigenvector selection from an alternative feedback design: There are also additional input arguments when we want to reproduce by projection (option 'p') eigenvector assignment performed by another gain (denoted fb0 obtained for example using µ- or LQR synthesis) in this case, we shall invoke fb_prop as follows:

  K = fb_prop(sys,0,pol_cl,'p',pol_fb0,fb0)

  in which pol_fb0 contains poles assigned by the (possibly dynamic) gain fb0.
• Multi-model case: Eigenstructure assignment constraints can be relative to more than one model. In that case, instead of fb_prop invoke defin_vw which has almost the same input arguments. This function computes assignment constraints instead of the feedback gain. For two models, the syntax will be

  [v1,w1,cv1] = defin_vw(sys1,pol_cl1,key1,def_pb1);
  [v2,w2,cv2] = defin_vw(sys2,pol_cl2,key2,def_pb2);
  K = [w1 w2]*inv([cv1 cv2]);

Example 2 : use of fb_prop with an observer

In this example, we consider the system of Example 1 from which two measurements are removed. Therefore, there are still 4 dominant modes but only 2 measurements. As all dominant modes must be treated, keeping in mind that modal control with proportional feedback is efficient when the number of measurements is at least equal to the number of treated modes, two measurements are missing. We shall use an observer of order 2 in order to replace the 2 missing measurements (note that we do not need a minimum order state observer that would be of order 4).

There are several functions available for observer design, see here. The function ob_ins optimizes the linear combinations of states that are observed in order to provide the best complement of the observation matrix (matrix usually denoted C). The function ob_gene is the dual counterpart of fb_prop. We shall use this function with the option 'p' which means projection of open-loop left eigenvectors. The benefits of such a strategy is discussed on page 51 of the manual.

So, a 2nd order observer is designed by projection of two open-loop left eigenvectors. The observer is connected to the system. So, 4 measurements (the 2 original measurements plus 2 observed signals) are available. Finally, 4 poles are shifted with "minimum energy". The design procedure reduces to:

   sys2 = [1 0 0 0;0 0 0 1]*rcamdata('lat',0,1);
   [u,t,pi] = ob_gene(sys2,[-1,-2],'p',[-0.18,-1.3]);
   sys3 = add_obs(sys2,u,t,1,pi);
   pol_ol = [-0.23+0.59j -0.18 -1.30];
   pol_cl = [-0.60+0.60j -0.60 -1.30];
   K3 = fb_prop(sys3,0,pol_cl,'m',pol_ol)

% Load a system with 6 states, 2 measurements and 2 inputs

   sys2 = [1 0 0 0;0 0 0 1]*rcamdata('lat',0,1);

% Compute open-loop eigenvalues

   eig(sys)

   ans =

     -1.3017
     -0.1837
     -0.2360 + 0.5954i
     -0.2360 - 0.5954i
     -6.6667
     -3.3333

% Observer designed such that two observed signals are
% added to the measurements: The observed signals are defined
% by projection of the left eigenvectors corresponding to the
% open-loop poles -0.18 and -1.30 and
%    -0.18  shifted to  -1
%    -1.30  shifted to  -2

   [u,t,pi] = ob_gene(sys2,[-1,-2],'p',[-0.18,-1.3]);

% This observer is connected to the system

   sys3 = add_obs(sys2,u,t,1,pi);

% sys3 considered as having 4 measurements with "C-matrix" [sys3.c;u]
%    -0.23+0.59j  replaced by  -0.60+0.60j
%    -0.18        replaced by  -0.60
%    -1.30        preserved

   pol_ol = [-0.23+0.59j -0.18 -1.30]; % This is def_pb
   pol_cl = [-0.60+0.60j -0.60 -1.30];

% Assignment by invoking fb_prop. Pole shifting is performed with
% minimum energy option ('m')

   K3 = fb_prop(sys3,0,pol_cl,'m',pol_ol)

   K3 =

      -1.8347    0.3168    1.3869    0.0241
      -2.0624   -0.4331    4.2717   -1.4181

% Results (eig_fb computes closed-loop eigenvalues)

   eig_fb(sys3,K3)

   ans =

     -6.7037
     -3.7334
     -2.0000
     -0.6000 + 0.6000i
     -0.6000 - 0.6000i
     -0.6000
     -1.0000
     -1.3000

% We recognize the observer poles (-1 and -2), the poles assigned
% by feedback (-0.6 + 0.6j, -0.6 and -1.3) and the non-dominant poles
% that have moved to -6.70 and -3.73.

Additional comments on observer-based control design.

• Transformation of a given dynamic controller to an equivalent observer-based controller: The function to be used is dfb2obs. The syntax is

  [u,t,pi,Q,K] = dfb2obs(sys,fb,0);

  in which fb is the dynamic feedback to be transformed. The output arguments u, t, pi have the same meaning as in Example 2. There is no constraint on the respective dynamic orders of sys and fb. The input argument 0 specifies that the slowest poles will be devoted to observation. Alternatively, -100 would mean that the poles the closest to -100 are assigned to observation. It is also possible to replace this argument by a vector which contains the closed-loop poles that must correspond to observation. The observer-based form of fb can be modified using a script similar to the one of Example 2, multi-model design can be used to improve robustness, order reduction becomes often straightforward and so on. In addition, the controller in this form can easily be initialized (initial value of observer state vector is u*x₀ where x₀ is the initial state).
• Direct observer-based feedback design: Two functions dfb_proj and dfb_ins perform design procedures that look like the one of Example 2. For example

  dfb = dfb_proj(sys,settling_time,damping_ratio,'p')
  

  However, it is very difficult to define a systematic modal control design procedure, therefore, these functions have some limitations. For example, dfb_proj does not work when the system to be controlled contains open-loop integrators.

Example 3 : elementary use of fb_tun

This function solves a Linear Quadratic Programming Problem. The input arguments of this function are two matrices:

• CSTR which defines linear constraints
• CRIT which defines a quadratic criterion.

In this example we push to the left the poles with real part larger than -0.5. For that purpose, an iterative procedure is written. At each step of this procedure, the poles with real part larger than -0.5 are detected by using the function sort_ev. Then, first order constraints corresponding to a minimum shift of 0.1 to the left, are defined using the function cstr_dp. Finally, the function fb_tun computes the first order gain variation corresponding to these constraints (LQ Programming problem in which the quadratic criterion is defined by default: CRIT = []). The correponding script is :

   sys = demodata(1);
   k0 = zeros(2,4);
   for ii = 1:100;
       lam = eig_fb(sys,k0);
       lamshift = sort_ev(lam,'r>-0.5&i>=0');
       if length(lamshift) == 0; disp('BREAK'), break,end;
       CSTR = cstr_dp([],sys,lamshift,'r','<',-0.1,k0);
       dk = fb_tun(CSTR,[]);
       k0 = k0+dk;
   end;

Let us examine more closely this script.

% The eigenvalues are shifted to the left until all real parts
% become less than -0.5

% Initial system and initial gain

   sys = demodata(1);
   k0 = zeros(2,4);

% Tuning loop

   for ii = 1:100;

% Eigenvalues to be shifted identified using sort_ev
% (real part > -0.5 and imaginary part positive or zero)

       lam = eig_fb(sys,k0);
       lamshift = sort_ev(lam,'r>-0.5&i>=0');
       if length(lamshift) == 0; disp('BREAK'), break,end;

% First order shifting constraint plugged into CSTR. Real part 
% of each pole in lamshift reduced of an amount equal to 0.1

       CSTR = cstr_dp([],sys,lamshift,'r','<',-0.1,k0);

% The gain variation dk is computed and the gain k0 is updated

       dk = fb_tun(CSTR,[]);
       k0 = k0+dk;

   end;

% Initial spectrum

   eig_fb(sys)

   ans =

     -0.1272 + 1.0734i
     -0.1272 - 1.0734i
     -0.9212
      0.0057
     -9.9000
    -10.1000          

% Results: final spectrum

   eig_fb(sys,k0)

   ans =

     -9.1834
     -9.8022
     -0.5835 + 0.5447i
     -0.5835 - 0.5447i
     -0.5087 + 0.1095i
     -0.5087 - 0.1095i

Additional comments on tuning.

Pole selection.

• The function sort_ev that is used to select poles to be treated at each step of the iterative tuning procedure in Example 3 can be used with complex logical expressions. For example if we want to detect the poles that satisfy one of the following conditions:
  • imaginary part between +1 and +2, real part less than -1 and damping ratio larger than 0.7,
  • imaginary part between +4 and +8 and real part between -1 and -2,
  we shall write:

  sort_ev(lambda,'(i>1&i<2&r<-1&d>0.7)|(i>4&i<8&r>-2&r<-1)')
  

• There are alternative functions for pole selection (see here). The function contr_ev can be used to eliminate eigenvalues with weak degree of controllability / observability. vicin_ev is an alternative tool which selects poles that belong to circles with given center and radius.

• The set of functions (see here): cstr_k (feedback gain constraints), cstr_eig (eigenstructure assignment), cstr_ev (eigenvalue assignment constraints) and cstr_qud (constraints for pole assignment in a domain) can be used like cstr_dp in the above example. This example illustrates tuning of the real part of some poles with cstr_dp. With the same function, it is also possible to tune the imaginary part and to constrain pole motion inside cones or on iso-damping lines.
• In order to avoid confusions with the functions for use with fb_dyn (see here), the functions relative to fb_tun (see here) have names of the form "cstr_.." or "crit_.." (those relative to fb_dyn have names of the form ".._cstr" or ".._crit").
• There are two kinds of satellite functions "cstr_.." : those which define absolute constraints (cstr_eig and cstr_k) and those which define first order constraints (cstr_dp, cstr_ev and cstr_qud). Some care must be exercised when constraints of both kinds are used together because absolute constraints might induce large gain variations that are in contradiction with first order approximations.
• The functions cstr_eig plus fb_tun permit the designer to perform non-iterative eigenstructure assignment (like using fb_prop, see Example 1), for example :

  CSTR = cstr_eig([],sys,pol_cl,'m',pol_ol,0,0);
  K1 = fb_tun(CSTR,[]);
  

  is equivalent to the invokation of fb_prop in Example 1. But here, cstr_k can be called before fb_tun for additional feedback gain constraints, multi-model design becomes also feasible by using several times cstr_eig with different models.

• The function fb_tun was written for multi-model control design : Each times one of the satellite functions "cstr_.." is invoked, the system considered might be different.
• There are some extensions to dynamic feedback gain tuning. The best approach for that consists of tranforming a given dynamic feedback to an observer-based feedback using dfb2obs (see here). In this case observer poles are not controllable, so cannot be tuned.

• For each first order constraint, the user must define the speed at which the feedback gain is tuned (for example, in Example 3, the aforementioned speed was the eigenvalue shift amount -0.1). If this speed is too large, the algorithm will not converge. If it is too small, a lot of iterations will be necessary. There is no rule for an a priori selection of the speed (usually, a few trials and errors are sufficient).
• Sophisticated optimization routines offer probably better convergence. The advantage of the proposed tuning procedure is its simplicity as it suffices to put together already existing tools in a very short script.

Example 4 : use of fb_dyn for multimodel "phase control"

Classical phase control consists of modifying the departure angles of some branches of the root locus. These directions are usually modified by tuning the phase of the system in given intervals of frequencies. On account of the use of the phase concept, this classical technique is obviously limited to SISO systems. In the manual on page 71, it is shown how the phase interpretation can be revisited in an eigenstructure assignment setting. This alternative interpretation leads to a generalization of "phase control" to MIMO and to multi-model control design.

The main function for "generalized phase control" is fb_dyn. This function solves a Linear Quadratic Programming Problem. The input arguments of this function are two matrices:

• CSTR which defines linear constraints
• CRIT which defines a quadratic criterion.

The structure of the dynamic feedback gain must be defined first by invoking the function str_cstr. These constraints initialize the matrix CSTR. Then, the function dp_cstr will be used in order to compute pole shifting constraints and to write them into the matrix CSTR. The feedback that is obtained must be multiplied by a scalar that must be tuned from zero to some positive value.

The illustrative example that is considered here consists of designing a feedback simultaneously relative to two systems. For both systems (each one has 2 pairs of complex conjugate poles), one complex pole must remain fixed and the other one must moves to the left along an horizontal tangent. Assuming that the two systems (sys1 and sys2) are in the work space, the design procedure reduces to the following script:

   CSTR = str_cstr(1,2,[-2+2*i;-3+3*i],1);
   CRIT = ktf_crit(CSTR,[0 0],i*[0 0.01 0.05 0.1 0.5 1]);

   pol1 = [ -0.1+3*i;-0.1+3.0*i;-1.0+2.0*i;-1.0+2.0*i];
   pol2 = [ -1.0+4*i;-1.0+4.0*i;-0.1+6.0*i;-0.1+6.0*i];
   dpol = [   -0.2  ;       0   ;     0   ;      0   ];
   CSTR = dp_cstr(CSTR,sys1,pol1,'riri','<===',dpol);
   CSTR = dp_cstr(CSTR,sys2,pol2,'riri','<===',dpol);

   fb = fb_dyn(CSTR,CRIT);

Root locus results:

Explanations relative to the above script.

% Two systems considered (1 input, 2 outputs).
% System 1 (poles -1+2i and -0.1+3i):
   den1 = [-1+2*i -1-2*i -0.1+3*i -0.1-3*i];
   sys1 = ss(zpk({[1];[2]},{den1;den1},[1;2]));
% System 2 (poles -1+4i and -0.1+6i): 
   den2 = [-1+4*i -1-4*i -0.1+6*i -0.1-6*i];
   sys2 = ss(zpk({[1];[3]},{den2;den2},[1;2]));

% Gain structure: size 1 x 2, two poles = -2+2*i and -3+3*i
% and denominator degree - numerator degree = 1
   CSTR = str_cstr(1,2,[-2+2*i;-3+3*i],1);

% Criterion: reference gain = [0 0] at listed frequencies
   CRIT = ktf_crit(CSTR,[0 0],i*[0 0.01 0.05 0.1 0.5 1]);

% Phase control constraints to be read column wise<
% pole -0.1+3i with real part variation < -0.2
% pole -0.1+3i with imaginary part variation = 0
% pole -1+2*i with real and imaginary part variations = 0
% Similar for the poles of sys2.
   pol1 = [-.1+3*i;-.1+3*i;-1+2*i ;-1+2*i ];
   pol2 = [ -1+4*i;-1+4*i ;-.1+6*i;-.1+6*i];
   fld1 = [ 'r'   ;  'i'  ;  'r'  ;   'i' ];
   fld2 = [ '<'   ;  '='  ;  '='  ;   '=' ];
   fld3 = [ -0.2  ;   0   ;   0   ;    0  ];
   [CSTR] = dp_cstr(CSTR,sys1,pol1,fld1,fld2,fld3);
   [CSTR] = dp_cstr(CSTR,sys2,pol2,fld1,fld2,fld3);

% Gain computation
   fb = fb_dyn(CSTR,CRIT);

% Root locus results
   figure; rlocus(sys1*fb,[0:-0.1:-2])
   figure; rlocus(sys2*fb,[0:-0.1:-2])

Example 5 : use of fb_dyn for multi-model eigenstructure assignment

Multi-model control laws cannot be designed in one shoot. It is recommended to design first a single-model nominal law, then to apply this law to other models (or use µ-analysis). Analysing other models permits the designer to identify one or more worst cases that must be treated together with the nominal system. When worst cases are treated, in order to avoid conflicting objectives, it is better to take into account the closed-loop properties of worst case models controlled by the initial feedback. For example, move the poles but try to preserve the assigned eigenvectors by using projection ideas.

The main function for multi-model eigenstructure assignment is fb_dyn. This function solves a Linear Quadratic Programming Problem. The input arguments of this function are two matrices:

• CSTR which defines linear constraints
• CRIT which defines a quadratic criterion.

The structure of the dynamic feedback gain must be defined first by invoking the function str_cstr. These constraints initialize the matrix CSTR. Then, the function eig_cstr will be invoked once for each selected worst case model. This function computes eigenstructure assignment constraints and write them into the matrix CSTR. The input arguments of eig_cstr are similar to those of fb_prop (without the second argument k0: see Example 1).

In the following illustrative example, it is assumed that the initial feedback fb0 has been designed relative to sys1 (assigning pol0) and, from analysis, the system sys2 has been identified as a worst case. Assuming that fb0, pol0, sys1 and sys2 are in the work space:

   CSTR = str_cstr(2,6,-8,0);
   CRIT = ktf_crit(CSTR,fb0,i*[0 0.01 0.05 0.1 0.5]);

   CSTR = eig_cstr(CSTR,sys1,pol0,'p',pol0,fb0);

   oldpol = [-2.37+4.40*i -1.14+3.0*i];
   newpol = [-4+4*i -3+3*i];
   CSTR = eig_cstr(CSTR,sys2,newpol,'p',oldpol,fb0);

   fb = fb_dyn(CSTR,CRIT);

Explanations relative to the above script (see Example 1 for the design of fb0).

% Generation of the nominal model 
   sys1 = rcamdata('lat',12,2,0.1);

% Initial proportional feedback computed for sys1
   pol0 = [-0.9+0.9i -0.8 -2 -1.5 -1];
   def_pb = [4 4 1 1 1];
   K0 = fb_prop(sys1,0,pol0,'n',def_pb);
   fb0 = ss([],[],[],K0);

% Analysis: the worst case model turns out to be sys2 below
   sys2 = rcamdata('lat',25,2,0.075);

% The system sys2 controlled by fb0 has two pairs of badly damped
% closed-loop poles which are -2.37 + 4.40i and -1.14 + 3.00i.
% These poles will be shifted to -4 + 4i and -3 + 3i. The
% initially assigned poles of sys1 will be re-assigned.

% Gain structure: size 2 x 6, a single pole at -8,
% and denominator degree - numerator degree = 0
   CSTR = str_cstr(2,6,-8,0);

% Criterion: reference gain fb0 at listed frequencies
   CRIT = ktf_crit(CSTR,fb0,i*[0 0.01 0.05 0.1 0.5]);

% For sys1 the initial assignment (by fb0) is preserved 
   CSTR = eig_cstr(CSTR,sys1,pol0,'p',pol0,fb0);

% For sys2 the poles assigned by fb0 are moved to the left
% and the corresponding eigenvectors are the projections of
% those assigned by fb0
   oldpol = [-2.37+4.40*i -1.14+3.0*i];
   newpol = [-4+4*i -3+3*i];
   CSTR = eig_cstr(CSTR,sys2,newpol,'p',oldpol,fb0);

% Dynamic feedback computation
   fb = fb_dyn(CSTR,CRIT);

% Analysis of the results (pol0 is in the first spectrum
% newpol in the second one):
    [eig_fb(sys1,fb) eig_fb(sys2,fb)]

    ans =

    -30.7701 +18.3849i -41.3772 +25.0282i
    -30.7701 -18.3849i -41.3772 -25.0282i
    -30.5619 +18.1003i -41.0257 +24.5727i
    -30.5619 -18.1003i -41.0257 -24.5727i
     -7.4785            -4.0000 + 4.0000i   <- assigned
     -6.5154            -4.0000 - 4.0000i
     -0.9000 + 0.9000i  -3.0000 + 3.0000i   <- assigned
     -0.9000 - 0.9000i  -3.0000 - 3.0000i
     -2.0475            -2.8696 + 0.2111i
     -2.0000            -2.8696 - 0.2111i
     -1.5647            -1.5690 + 1.0900i
     -1.5000            -1.5690 - 1.0900i
     -1.0000            -0.5829
     -0.8000            -0.3766

Example 6 : use of fb_dyn for controller order reduction

The proposed method for controller order reduction consists of:

• Identifying the dominant eigenstructure assigned by a given high order controller.
• Designing a new controller that only assigns this dominant eigenstructure.

The main function for controller order reduction (i.e. re-design) is fb_dyn. This function solves a Linear Quadratic Programming Problem. The input arguments of this function are two matrices:

• CSTR which defines linear constraints
• CRIT which defines a quadratic criterion.

The structure of the dynamic feedback gain must be defined first by invoking the function str_cstr. These constraints initialize the matrix CSTR. Then, the function eig_cstr is invoked in order to define constraints corresponding to re-assigning the dominant eigenstructure. This function computes these constraints and write them into the matrix CSTR. The input arguments of eig_cstr are similar to those of fb_prop (without the second argument k0: see Example 1).

The illustrative example considered below was proposed by Caroline Chiappa in [Chiappa et al, Improvement of the robustness of an aircraft autopilot designed by a µ-synthesis technique, CESA'98].

The initial controller was designed using a µ-synthesis approach detailed in [Bennani and Looye, Chapter 22 of Robust flight control: a design challenge, LNCIS 224]. The considered system is the longitudinal motion of the RCAM (the linearized models used here have more states and measurements than the models available by invoking the RMCT function rcamdata, therefore, augmented models are provided in a .zip file). The controller architecture is depicted in the following figure. For notation, see the aforementioned papers.

The initial feedback fb0 (K(s) in the figure) has been designed relative to the system sys0 (RCAM in the figure). Assuming that fb0, sys0 are in the work space:

   sysbf = feedback(sys0,fb0,1);
   p_dom_bf = plot_res(sysbf,1,1,inf);
   p_dom = p_dom_bf(1:5);

   beta=[-10 -0.015];
   [CSTR] = str_cstr(3,6,beta);

   omega=[.01:.01:0.1]*sqrt(-1);
   CRIT = ktf_crit(CSTR,fb0,omega);

   CSTR = eig_cstr(CSTR,sys0,p_dom,'p',p_dom,fb0);

   fb = fb_dyn(CSTR,CRIT,0.0001);

Explanations relative to the above script.

% The system sys0, and the initial feedback fb0 are loaded
   load ex6_data

% Analysis for identifying dominant modes
   sysbf = feedback(sys0,fb0,1);
   [p_dom_bf,res] = plot_res(sysbf,1,1,inf);
   format short e
   [p_dom_bf,res]

    ans =

      -6.4634e-01 + 2.9574e-01i -4.3401e+01
      -1.1819e+00                3.2808e+01
      -2.5150e-01 + 6.2944e-01i  1.2327e+01
      -8.3603e-01               -2.8614e+00
      -4.0363e-01 + 6.1812e-01i  1.0316e+00
      -6.5799e+00 + 5.8447e+00i  1.3467e-01
      -7.1421e-01               -3.0885e-02
      -3.1618e-01 + 3.1629e-01i  1.8490e-03
      -1.5050e-01                1.0791e-03
      -1.4950e-01               -1.0737e-03
      -1.1248e-01 + 3.1673e-02i -4.3618e-05
      -1.4144e-01 + 1.4136e-01i  2.7219e-06
      -1.1143e+02                2.2292e-06
      -1.1042e+01                2.3912e-07
      -1.0940e+02 + 1.1127e+02i -1.3514e-07
      -1.2821e+02                9.1333e-10
      -1.9876e+03                1.9283e-15
      -6.7000e-01                5.6541e-22  

% For example, consider 5 poles (residuals >1, right column)
   p_dom = p_dom_bf(1:5);

% Feedback gain structure definition: size 3 x 6, two poles
% -100 and -0.015. One of these values must be close to the
% origin in order to preserve the integral effect of fb0
   beta=[-10 -0.015];
   [CSTR] = str_cstr(3,6,beta);

% Criterion: reference gain fb0 at listed frequencies
   omega=[.01:.01:0.1]*sqrt(-1);
   CRIT = ktf_crit(CSTR,fb0,omega);

% Constraints for re-assignment of dominant modes (poles
% re-assigned here at the same location because the 3rd and 
% 5th arguments are equal) 
   CSTR = eig_cstr(CSTR,sys0,p_dom,'p',p_dom,fb0);

% Dynamic feedback computation
   fb = fb_dyn(CSTR,CRIT,0.0001);

% Step responses analysis
   figure
   sys_k = fb;
   plotstep
   figure
   sys_k = fb0;
   plotstep

Additional comments on order reduction.

Selection of closed-loop dominant modes

• The above example is somewhat misleading because dominant modes selection turns out to be extremely easy as it suffices to use only once plot_res between any input-output pair. The function lsim_mod gives more information and might be used preferably to plot_res.
• Using lsim_mod some closed-loop modes appearing to be dominant in the modal simulations might cancel each other, therefore are not dominant. In fact dominant modes selection cannot be fully automatized as it requires generally expertise from the designer.
• In the open-loop case, dominant mode identification is almost always obvious. So, closed-loop dominant poles can often be partly identified systematically by a "root locus" corresponding to a feedback ρfb0 with ρ varying form 0 to 1(fb0 denoting the feedback gain to be reduced): If there is no branch crossing, branches starting from dominant open-loop poles generally lead to final points that are dominant closed-loop poles. Note that the number of dominant modes is often larger in closed-loop than in open-loop.

Controller poles

• The choice of the controller poles does not affect significanlty the quality of the results provided that a few basic rules are followed. Briefly, denominator poles must have dynamics corresponding to the dominant modes that are re-assigned (e.g., do not consider all denominator poles around -0.1 if one of the modes to be re-assigned is as "fast" as -10. If there is an integral effect, consider a "slow" pole.).
• In most systems we have D=CB=0 (standard notation). In this case, it can be shown that the sum of the controller poles and of the open-loop system poles is constant. It means that pushing some open-loop poles to the left will result in controller poles (or other non-dominant poles) to be shifted of the same amount to the right. Clearly, this rule helps to select denominator poles.
• Denominator poles have limited effect because in most cases, using dynamic controllers corresponds to using derivatives of measurments as additional measurements. In such cases, the denominators make derivatives realizable as pseudo-derivatives. The corresponding bandwidth must depends on disturbances (and expected roll-off properties).
• In the MIMO case, in order to optimize controller order reduction, it is better to use the same denominator for all the controller entries. Explanation: each pole will be "repeated" at least once but if it appears several times in the same row (or column) it will be "repeated" only once, see minimal realization theory.

Frequency domain constraints

• Example 6 shows how to preserve time domain specifications. If there are also frequency domain properties (like a stop band filter) to be recovered, it is suggested to use the function add_cstr. For example if we want the distance between fb0 and the gain being computed to satisfy:
  • distance of entry (1,2) less than 0.1 at 1.1 rd/s
  • distance of entry (3,2) equal to 0.0 at 2.2 Rd/s
  we shall insert in a script before fb_dyn is invoked

  CSTR = add_cstr(CSTR,[1 2;3 2],[1.1;2.2],['<';'='],[0.1;0.0],fb0)

  It is also possible to use constraints that do not depend on fb0, for example if we want all entries to be less than 0.01 at 20 rd/s, we shall use

  CSTR = add_cstr(CSTR,'all',20,'<',0.01)

  Note that at 0 Rd/s (and only at this frequency) we can define a constraint of the form "larger than". For example, entry (2,2) with DC gain larger than 10:

  CSTR = add_cstr(CSTR,[2 2],0,'>',10)