Why does augw(P, W1, W2, W3) give a much lower order generalized plant when compared to the by hand construction?

I have found that constructing a generalized plant with the built in function augw(P, W1, W2, W3) leads to a much lower order plant than just constructing one by hand as specified in its reference page:

https://www.mathworks.com/help/robust/ref/dynamicsystem.augw.html#mw_4dda95b9-fe05-4a4d-8930-9e60892694e1

This is the case even when considering the minimal realization of the by hand construction. Just to demonstrate this, I generated a bunch of random test models and weights in the TestModels.mat file linked above, and you can see the discrepancy using the script below:

load('TestModels.mat');
TestModel=TestStruct.TestModel; %Test plant
TestW1=TestStruct.TestW1;       %Test W1 weight
TestW2=TestStruct.TestW2;       %Test W2 weight
TestW3=TestStruct.TestW3;       %Test W3 weight
Ze=zeros(2,2);

%Formulating the generalized plant exactly as described in augw() reference
GgenByHand=[TestW1       -TestW1*TestModel; 
            Ze            TestW2;
            Ze            TestW3*TestModel;
            eye(2,2)     -TestModel];
 
 
 GgenAugW=augw(TestModel,TestW1,TestW2,TestW3); %Formulating the generalized plant using augw() function
 
 %GgenAugW has a much smaller order even when considering the minimal realization of the by hand generalized plant
 size(GgenAugW) 

State-space model with 8 outputs, 4 inputs, and 23 states.

 size(minreal(GgenByHand))

19 states removed.
State-space model with 8 outputs, 4 inputs, and 42 states.

I was wondering, why is this the case? Again going off of the augw reference, it seems like the by hand construction is exactly what it does, but I don't see how they could have a lower order system when compared to even the minimal realization of that generalized plant. Doing some testing, it seems like a lot of the time the minimal realization order and the augw() order both agree, but unfortunately for my actual plant model this is not the case.

[Vinh]

Oh thank you Walter Roberson for the helpful edits.

Answer 1

[Paul]

Hi Vinh,

The 'by-hand' construction is nearly certainly not what augw actually does. You can look at the augw code; maybe it will be clear how it works.

load('TestModels.mat');
TestModel=TestStruct.TestModel; %Test plant
TestW1=TestStruct.TestW1;       %Test W1 weight
TestW2=TestStruct.TestW2;       %Test W2 weight
TestW2 = ss(TestW2);
TestW3=TestStruct.TestW3;       %Test W3 weight
Ze=zeros(2,2);
%Formulating the generalized plant exactly as described in augw() reference
GgenByHand=[TestW1       -TestW1*TestModel; 
            Ze            TestW2;
            Ze            TestW3*TestModel;
            eye(2,2)     -TestModel];
 
GgenAugW=augw(TestModel,TestW1,TestW2,TestW3); %Formulating the generalized plant using augw() function
 

As noted these models have a different number of states

order(GgenByHand)

ans = 61

order(GgenAugW)

ans = 23

The former has as many states as is needed to construct each element of that matrix

order(TestW1) + order(TestW1) + order(TestModel) + order(TestW2) + order(TestW3) + order(TestModel) + order(TestModel)

ans = 61

while the latter has only as many states needed to construct the augmented system

order(TestModel) + order(TestW1) + order(TestW2) + order(TestW3)

ans = 23

In theory, minreal should (I think) be able to get a realization with the same number of sates as in GgenAugW, but in practice it might not be so easy.

As you've seen, minreal can get closer with the default tolerance

order(minreal(GgenByHand))

19 states removed.

ans = 42

Play with the tolerance to see if we can get closer. A tolerance of 1e-6 is pretty large

sysmin = minreal(GgenByHand,1e-6); order(sysmin)

23 states removed.

ans = 38

but doesn't get all the way there. At least the frequency response of sysmin still matches.

bode(GgenAugW,GgenByHand);

If you want to construct the model by hand, it's better to use interconnection functions. Personally, I prefer connect, but other approaches are possible. For example

set(TestModel,'InputName','u','OutputName','y');
set(TestW1,   'InputName','e','OutputName','z1');
set(TestW2,   'InputName','u','OutputName','z2');
set(TestW3,   'InputName','y','OutputName','z3');

eblk = sumblk('e = r - y',2);

sys = connect(TestModel,TestW1,TestW2,TestW3,eblk, ...
    [eblk.InputName(1:2); TestModel.InputName], ...
    [TestW1.OutputName; TestW2.OutputName; TestW3.OutputName; eblk.OutputName]);
size(sys)

State-space model with 8 outputs, 4 inputs, and 23 states.

size(GgenAugW)

State-space model with 8 outputs, 4 inputs, and 23 states.

Compare

figure
bode(GgenAugW,sys)

Here's another approach from the perspective of input/output relationshships

%[w;u;y]      = [I 0;0 I;0 G] * [w;u]
%[e;u;y]      = [I 0 -I;0 I 0;0 0 I] * [w;u;y]
%[z1;z2;z3;e] = [W1 0 0;0 W2 0; 0 0 W3;I 0 0] * [e;u;y]
I2 = ss(eye(2)); z2 = ss(zeros(2));
sys = [append(TestW1,TestW2,TestW3);[I2 z2 z2]] * [I2 z2 -I2;z2 I2 z2;z2 z2 I2] * [I2 z2;z2 I2;z2 TestModel];
order(sys)

ans = 23

figure
bode(GgenAugW,sys)

Check the doc pages on Model Interconnection, particularly Preventing State Duplication in System Interconnections

[Vinh]

Thanks Paul!

[Paul]

You're quite welcome. I edited the answer to add another by-hand approach.