modred() Scaled Output

Hi MathWorks Community,

I am a Flight Controls Engineer and often use the modred() function to reduce a higher-order aircraft state-space model to a 6-DOF model for simplified design and analysis. I recently had an issue where modred() operations in newer versions of MATLAB (e.g. R2023b) no longer gave the same result as operations in older versions (e.g. R2018b).

This was a problem that for a long time I was unable to debug until my colleague recently noticed that the MathWorks documentation for newer MATLAB versions (e.g. R2023b) states "modred returns a scaled version of this realization. To disable this scaling, set sys.Scaled to true before eliminating the states."

This fixed the issue, but I wanted to ask why the default output would be a "scaled version of the realization" and not the result at the original scale. I can see that the scaling preserves the input-to-output relationship and provides robustness to perturbation of the parameters (https://www.mathworks.com/help/control/ug/scaling-models-to-maximize-accuracy.html), but in many cases (especially for physical systems) it is useful to preserve the original scaling by default such that the state-space derivatives are valid physical parameters corresponding to the original states, inputs, and outputs.

Thanks,

Sam

[Paul]

Hi Sam,

Can you upload a .mat file with an example system (use the paper clip icon on the Insert menu) and post the code that illustrates the problem?

[Samuel Nadell]

Hi @Paul,

Please find the example now attached. Let me know if that helps.

Thanks,

Sam

Answer 1

[Paul]

Hi Sam,

I can't explain why TMW chose to change the behavior of modred.

Having said that, I have some observations about your example.

dbtype Example_20240221_FINAL.m

1     % Example for modred() MathWorks Answers post
2     % 
3     % Refs:
4     % https://www.mathworks.com/help/control/ref/statespacemodel.modred.html
5     % * 'MatchDC' (default): Enforce matching DC gains. The state-space matrices are recomputed as described in Algorithms.
6     % 
7     % * elim can be a vector of indices or a logical vector commensurate with X where true values mark states to be discarded. 
8     %   This function is usually used in conjunction with balreal. Use balreal to first isolate states with negligible 
9     %   contribution to the I/O response. If sys has been balanced with balreal and the vector g of Hankel singular values has 
10    %   M small entries, you can use modred to eliminate the corresponding M states. For example:
11    % 
12    %       [sys,g] = balreal(sys)  % Compute balanced realization
13    %       elim = (g<1e-8)         % Small entries of g are negligible states
14    %       rsys = modred(sys,elim) % Remove negligible states
15    
16    clc;
17    clear;
18    close all;
19    
20    
21    %% Formulate the reduced-order systems
22    % Use XV-15 Hover from Tischler 2017
23    % Inputs: aileron [deg], rudder [deg]
24    % States: v [ft/s], p [deg/s], r [deg/s], phi [deg]
25    % Outputs: v [ft/s], p [deg/s], r [deg/s], phi [deg], ay [ft/sec^2]
26    load('Ex3_XV15_Hover_Model.mat')
27    sys = ss(A,B,C,D);
28    
29    % Check which states have small enough impact to eliminate
30    [sys_bal,g_bal] = balreal(sys);
31    
32    % Modred 1 (scale automatically modified by MATLAB by default)
33    sys_red = modred(sys,4);
34    
35    % Modred 2 (back to original scale)
36    sys.Scaled = true;
37    sys_red_scaled = modred(sys,4);
38    
39    % Bode
40    bodeplot(sys,'k-',sys_red,'r--',sys_red_scaled,'g.')
41    legend('full-order sys','modred (default)','modred (orig scale)')
42    
43    
44    %% Test impact on control allocation formulation
45    % m (desired moments: pdot, rdot) = CA * u (control inputs: ail, rud)
46    % --> u = inv(CA) * m
47    
48    % Original B-matrix
49    disp('The full-order system B-matrix has off-axis terms, indicating coupling:')
50    sys.B
51    
52    % Original scale
53    B_red_scaled = sys_red_scaled.B(2:3,:);
54    sys_scaled_decoupled = sys * inv(B_red_scaled);
55    disp('Using the .Scaled = true B-matrix in the control allocation properly decouples the')
56    disp('control inputs and has the proper scale (since it contains the identity matrix):')
57    sys_scaled_decoupled.B
58    % sys_scaled_decoupled.B =
59    %         u1      u2
60    % x1  0.7833       0
61    % x2       1       0
62    % x3       0       1
63    % x4       0       0
64    % 
65    % Since rows 2 and 3 are the identity matrix, it is clear that the controls
66    % have been successfully decoupled in the full-order system and are the 
67    % proper scale!
68    
69    % However, if do not set .Scaled = true, the combined system is decoupled
70    % but does not have the proper scaling
71    B_red = sys_red.B(2:3,:);
72    sys_decoupled = sys * inv(B_red);
73    disp('Using the default modred B-matrix results in a control allocation')
74    disp('not at the proper / physical scale (since the values are not ones):')
75    sys_decoupled.B
76    % sys_decoupled.B =
77    %          u1       u2
78    % x1  0.01224        0
79    % x2  0.01562        0
80    % x3        0     0.25
81    % x4        0        0
82    
83    

I see that the code is making a call to balreal, but the results from balreal aren't used.

Load in the model, create the ss object, and assign state, input, and output names in accordance with the comments in the .m file.

%Example_20240221_FINAL

load Ex3_XV15_Hover_Model
sys = ss(A,B,C,D,'StateName',{'v','p','r','phi'},'InputName',{'da','dr'},'OutputName',{'v','p','r','phi','ay'});

Something doesn't look right about sys. According to the comments, the fourth state is phi and the second state is p. The A matrix is correct in that it shows phidot = p. But it looks like the names of the fourth and fifth output variable are reversed. Clearly the fifth output is phi, not ay.

sys

sys =
 
  A = 
                v          p          r        phi
   v     -0.09785     -1.503          0      32.17
   p    -0.004377    -0.2362          0          0
   r     0.000716     0.0387    -0.1416          0
   phi          0          1          0          0
 
  B = 
              da        dr
   v    -0.04524         0
   p    -0.05776         0
   r    0.005908   0.01187
   phi         0         0
 
  C = 
               v         p         r       phi
   v           1         0         0         0
   p           0         1         0         0
   r           0         0         1         0
   phi  -0.09785    -1.503         0  0.004002
   ay          0         0         0         1
 
  D = 
              da        dr
   v           0         0
   p           0         0
   r           0         0
   phi  -0.04524         0
   ay          0         0
 
Continuous-time state-space model.

Call modred to eliminate the fourth state, which is also the fifth output. I suppose a state that is an output could be eliminated, but does it make sense to do so?

sys_red = modred(sys,4);

Now, even though the fourth state was supposed to be eliminated, the resulting output of modred still has four states:

size(sys_red)

State-space model with 5 outputs, 2 inputs, and 4 states.

so that fourth state was not eliminated. But something happened because sys_red is returned in descriptor form

sys_red.E

ans = 4×4

1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0

Now scale the system and try modred

sys.Scaled = true;
sys_red_scaled = modred(sys,4);

Again, no state was actually eliminated.

size(sys_red_scaled)

State-space model with 5 outputs, 2 inputs, and 4 states.

In fact, the A/B/C/D matrices of the output are just that of the input, which is why your control allocation scheme appears to work.

[sys_red_scaled.A sys_red_scaled.B;sys_red_scaled.C sys_red_scaled.D] - [sys.A sys.B;sys.C sys.D]

ans = 9×6

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

But, for some reason, sys_red_scaled is in descriptor form

sys_red_scaled.E

ans = 4×4

1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0

which I'll guess is the reason for the difference in this Bode plot (keeping in mind that sys is now scaled)

bodeplot(sys(1,1),sys_red_scaled(1,1))

At this point, I'm not sure what modred is doing for this example, but the results sure look odd (pending further investigation, I'm intrigued).

p.s. If we continue this thread, please don't post your code as a .m file attachment. Just copy/paste it directly into the comment and make sure to use the code formatting for the actual code (select the code and then click the leftmost icon in the Code ribbon).

[Paul]

Additional thoughts on this subject.

Load the model, assign the names, assume that ay is the fourth ouput an phi is the fifth

load Ex3_XV15_Hover_Model
sys = ss(A,B,C,D,'StateName',{'v','p','r','phi'},'InputName',{'da','dr'},'OutputName',{'v','p','r','ay','phi'})

sys =
 
  A = 
                v          p          r        phi
   v     -0.09785     -1.503          0      32.17
   p    -0.004377    -0.2362          0          0
   r     0.000716     0.0387    -0.1416          0
   phi          0          1          0          0
 
  B = 
              da        dr
   v    -0.04524         0
   p    -0.05776         0
   r    0.005908   0.01187
   phi         0         0
 
  C = 
               v         p         r       phi
   v           1         0         0         0
   p           0         1         0         0
   r           0         0         1         0
   ay   -0.09785    -1.503         0  0.004002
   phi         0         0         0         1
 
  D = 
              da        dr
   v           0         0
   p           0         0
   r           0         0
   ay   -0.04524         0
   phi         0         0
 
Continuous-time state-space model.

Looking at the D matrix, is it correct that the aileron has a direct feed to ay but dr does not?

Anywya, just considering the math....

Call modred incorrectly (don't include the second argument)

try
    sys_red1 = modred(sys)
catch ME
    ME.message
end

ans = 'The "xelim" command requires two or more input arguments. Type "help xelim" for more information.'

It looks like in 2023b modred (which is no longer recommended according to the doc page) just ships over to xelim.

Run modred on the system, as above the output is in descriptor form

sys_red1 = modred(sys,4)

sys_red1 =
 
  A = 
              v          p          r          ?
   v   -0.09785  -0.005873          0     0.2514
   p     -1.121    -0.2362          0          0
   r    0.01146   0.002418    -0.1416          0
   ?          0        0.5          0          0
 
  B = 
            da        dr
   v  -0.01131         0
   p    -3.697         0
   r   0.02363   0.04748
   ?         0         0
 
  C = 
                v          p          r          ?
   v            4          0          0          0
   p            0    0.01562          0          0
   r            0          0       0.25          0
   ay     -0.3914   -0.02349          0  0.0001251
   phi          0          0          0    0.03125
 
  D = 
              da        dr
   v           0         0
   p           0         0
   r           0         0
   ay   -0.04524         0
   phi         0         0
 
  E = 
      v  p  r  ?
   v  1  0  0  0
   p  0  1  0  0
   r  0  0  1  0
   ?  0  0  0  0
 
Continuous-time state-space model.

If you check the Algorithms section of modred for the MatchDC case, we see that the A22 of the partitioned system must be invertible. In our case A22 = A(4,4) = 0, and I suspect that's the root of the behavior we're seeing. Make A(4,4) a small value and try again.

sys.A(4,4) = -1e-5;
sys_red2 = modred(sys,4)

sys_red2 =
 
  A = 
              v          p          r
   v   -0.09785  1.257e+04          0
   p     -1.121    -0.2362          0
   r    0.01146   0.002418    -0.1416
 
  B = 
            da        dr
   v  -0.01131         0
   p    -3.697         0
   r   0.02363   0.04748
 
  C = 
              v        p        r
   v          4        0        0
   p          0  0.01562        0
   r          0        0     0.25
   ay   -0.3914    6.229        0
   phi        0     1562        0
 
  D = 
              da        dr
   v           0         0
   p           0         0
   r           0         0
   ay   -0.04524         0
   phi         0         0
 
Continuous-time state-space model.

Now phi is actually eliminated. Maybe modred, or actually xelim, should have at least put out a warning for the original case, and/or the relevant doc pages should have a discussion of what happens for MatchDC when A22 is not invertible.

Whether or not sys_red2 is an accurate enough representation of sys is a different question.

I believe the Algorithm sections of the doc pages for modred and xelim have a typo for the MatchDC case. It seems to me the B-matrix for the reduced system should be:

B1 - inv(A22)*B2 as opposed to: B1 - A12*inv(A22)*B2

[Samuel Nadell]

Hi Paul,

Thanks for the detailed observations and sorry it probably wasn't the best example. (By the way, to answer your questions, yes ay and phi should be switched in the outputs and it is correct that in the D-matrix Ydr = 0 (zero ay response to rudder in this model)). I guess to summarize the answer to my primary question, we "can't explain why TMW chose to change the behavior of modred" and users should make sure to set "sys.Scaled = true" in newer versions of MATLAB if they want to preserve the original scaling such that the state-space derivatives are valid physical parameters corresponding to the original states, inputs, and outputs. Nevertheless, your observations were insightful on use of modred and I'll close this question for now.

Best regards,

Sam

[Paul]

I'm goint to retract my statement that A22 must be invertible. Back in 2020b, the modred doc had a section on Limitations that stated that A22 had to be invertible for the MatchDC case (continuous time). But that limitation was removed in 2021a, so it's reasonable to assume that the function was changed in 2021a to allow for A22 to be singular, though I didn't see anything to this effect in the 2021a release notes. Too bad the current documentation doesn't explain exactly what happens when A22 is singular. Actually, the criterion must be something like close to singular because in the example we're working with this in this question if we change A(4,4) to a small, nonzero number

sys.A(4,4) = -1e-10;

I found that the result from modred came back in descriptor form.