Step() not working as expected with USS state space (Uncertain Sys)

I'm trying to figure out what step() is doing with uss (uncertain) discretized state space equations, for 2nd-order systems.

The response seems to sometimes blow up exponentially, or sometimes have weird transient behavior inconsistent with the dynamic. Any open-loop 2nd-O system is inherently stable and state should not grow unboundedly, or have odd dynamic response.

My process:

1) create uncertain (uss) state space from continuous version

2) plot continuous version, and discrete version

The discrete versions seem to be doing something odd: they have long-term transients inconsistent with the initial dynanics response. The cts versions don't have this (they behave as expected).

What's occurring here?

Ts = 1/50e3;

% Any values are a stable system with expected 2nd-O response
P = {};
P.J_kgPm2 = 1.0000e-05;
P.b = 3.2000e-03;
P.k_kgPm = 2.5266e+00;

multJ_u = ureal('multJ_u', 1, 'percent', 10);
multB_u = ureal('multB_u', 1, 'percent', 5);
multK_u = ureal('multK_u', 1, 'percent', 5);
P_u = {};
P_u.J_kgPm2 = P.J_kgPm2 * multJ_u;
P_u.b = P.b * multB_u;
P_u.k_kgPm = P.k_kgPm * multK_u;

% Make continuous-time uss system
A = [0, 1; -P_u.k_kgPm / P_u.J_kgPm2, -P_u.b / P_u.J_kgPm2];
B = [0; 1 / P_u.J_kgPm2];
C = [1 0];
D = [0];
sys = ss(A, B, C, D);
sys = ss(A, B/dcgain(sys), C, D, ...
    'StateName', {'theta', 'angVel'}, ...
    'InputName', {'u'}, ...
    'OutputName', {'theta'});

% Make discrete-time uss system
% c2d() doesn't support uss; make it manually from eg 1000 samples
numModels = 1000;
sysd_u = c2d(usample(sys, numModels), Ts, 'tustin');
sysd_u_nom = c2d(sys.NominalValue, Ts, 'tustin');
% Recreate uss-type by recombining disc models. 
% Note: 0 is to keep same # of states
sysd = ucover(sysd_u, sysd_u_nom, 0);
% Do it for long time, to see divergence
figure; 
step(sys, 0.3); grid on;
title('Cts. Expected steady-state behavior.')

figure;
step(sysd, 3); grid on;
title('Disc. Note the transients.')

figure;
step(sysd, 0.2); grid on;
title('Disc Zoom. Note the inconsistent transient behavior.')

For a case that blows up, here's an external-disturbance-augmented plant system. For a reference response, Cts version is fine, but Disc version blow up. Again, this is open loop, so responses should always settle.

% External dist force plant augmentation
A_aug = [ sys.A , sys.B;
    zeros(1, 2), 1];
B_aug = [sys.B; 0];
C_aug = [sys.C, 0];
% Input on the f-level
F_aug = [0; 0; 1 / P_u.J_kgPm2];
sysAugF = ss(A_aug, [B_aug F_aug], C_aug, 0);

% Make disc version
spread = c2d(usample(sysAugF, numModels), Ts, 'tustin');
nom = c2d(sysAugF.NominalValue, Ts, 'tustin');
sysAugFd = ucover(spread, nom, 0);

% Only ref response
figure;
subplot(2,1,1);
step(sysAugF(1), 0.05); title('Cts Aug Plant')
subplot(2,1,2);
step(sysAugFd(1), 0.05); title('Disc Aug Plant')
sgtitle('Augmented plant')

Answer 1

[Paul]

Hi John,

Disclaimer: I'm not really a user of the Robust Control Toolbox.

I'm not quite sure how to interpret the discretization of a uss model. Will have to think about that some more.

I think the fundamental issue is that the zero'th order model in ucover is too conservative.

Ts = 1/50e3;
% Any values are a stable system with expected 2nd-O response
%P = {};
P.J_kgPm2 = 1.0000e-05;
P.b = 3.2000e-03;
P.k_kgPm = 2.5266e+00;
multJ_u = ureal('multJ_u', 1, 'percent', 10);
multB_u = ureal('multB_u', 1, 'percent', 5);
multK_u = ureal('multK_u', 1, 'percent', 5);
%P_u = {};
P_u.J_kgPm2 = P.J_kgPm2 * multJ_u;
P_u.b = P.b * multB_u;
P_u.k_kgPm = P.k_kgPm * multK_u;
% Make continuous-time uss system
A = [0, 1; -P_u.k_kgPm / P_u.J_kgPm2, -P_u.b / P_u.J_kgPm2];
B = [0; 1 / P_u.J_kgPm2];
C = [1 0];
D = [0];
sys = ss(A, B, C, D);
sys = ss(A, B/dcgain(sys), C, D, ...
    'StateName', {'theta', 'angVel'}, ...
    'InputName', {'u'}, ...
    'OutputName', {'theta'});
numModels = 100; % reduced to meet Answers runtime requirement
sysd_u = c2d(usample(sys, numModels), Ts, 'tustin');
sysd_u_nom = c2d(sys.NominalValue, Ts, 'tustin');
% Recreate uss-type by recombining disc models. 
% Note: 0 is to keep same # of states
[sysd,Info] = ucover(sysd_u, sysd_u_nom, 0);

Uncertain state-space model with 1 outputs, 1 inputs, 2 states, and 1 blocks.

relerr = (sysd_u_nom-sysd_u)/sysd_u_nom;
figure
bodemag(relerr,Info.W1,'r')

We see that the constant weighting function allows for considerably more uncertainty at low frequencies.

Instead, try a second order uncertainty (I saw the comment about not wanting to increase the number of states, but was curious to try anyway)

[sysd,Info] = ucover(sysd_u, sysd_u_nom, 2);

Uncertain state-space model with 1 outputs, 1 inputs, 4 states, and 1 blocks.

relerr = (sysd_u_nom-sysd_u)/sysd_u_nom;
figure
bodemag(relerr,Info.W1,'r')

Now the uncertainty model hugs the upper bound of the error, and the responses look better, but it looks like there may still be some cases with undesirable transients. I suspect this is the price to pay of using the frequency dependent uncertainty model that probably allows for uncertainty that is not realizable by the uncertainty in the physical parameters.

figure;
step(sysd, 3); grid on;

figure;
step(sysd, 0.2); grid on;

[John]

I see. Thanks @Paul.

"I suspect this is the price to pay of using the freuqency dependent uncertainty model that probably allows for uncertainty that is not realizable by the uncertainty in the physical parameters."

Hmmm, agreed; that would be consistent...

I'm not quite sure how to proceed, then, because the system needs to maintain two states, so that I can create discretized controllers with known states (to expand into my other known controller structures).

Does anyone know if there's another way to create a discretized Uncertain model aside from ucover(), since c2d() doesn't seem suitable or uss?

[Paul]

I'll offer this up for you to think about. It looks promising, but more analysis required IMO. Buyer beware.

Define the system

Ts = 1/50e3;
% Any values are a stable system with expected 2nd-O response
%P = {};
P.J_kgPm2 = 1.0000e-05;
P.b = 3.2000e-03;
P.k_kgPm = 2.5266e+00;
multJ_u = ureal('multJ_u', 1, 'percent', 10);
multB_u = ureal('multB_u', 1, 'percent', 5);
multK_u = ureal('multK_u', 1, 'percent', 5);
%P_u = {};
P_u.J_kgPm2 = P.J_kgPm2 * multJ_u;
P_u.b = P.b * multB_u;
P_u.k_kgPm = P.k_kgPm * multK_u;
% Make continuous-time uss system
A = [0, 1; -P_u.k_kgPm / P_u.J_kgPm2, -P_u.b / P_u.J_kgPm2];
B = [0; 1 / P_u.J_kgPm2];
C = [1 0];
D = [0];
sys = ss(A, B, C, D);
sys = ss(A, B/dcgain(sys), C, D, ...
    'StateName', {'theta', 'angVel'}, ...
    'InputName', {'u'}, ...
    'OutputName', {'theta'});

Because we are using the Tustin rule for the transformation, we can just hard code the transformation equations on the state matrices of the continuous time uss model

sysd = tempbilinear(sys,Ts);

Compare step and frequency responses.

figure
step(sys,sysd)

figure
bode(sys,sysd)

I'm hesitant about this approach because, at least in part, it's not supported by c2d and maybe there's a good reason for that.

I suppose that one way to test this would be to realize sys with a particular values of the uncertainties as an ss, then do c2d on that ss, then compare that result to the realization of sysd as an ss by substituting into sysd the identical values of uncertainties that went into sys. I assume there's a way to do that. As I said previously, I'm not that familiar with RCT functionality.

Here's the transformation function, based on the algorithm used in bilinear.

function sysd = tempbilinear(sysc,Ts)

A = sysc.A;
B = sysc.B;
C = sysc.C;
D = sysc.D;

lambda = 1/Ts;
Ad = (eye(size(A)) - A/2/lambda) \ (eye(size(A)) + A/2/lambda);
Bd = 1/sqrt(lambda) * ((eye(size(A)) - A/2/lambda) \ B);
Cd = 1/sqrt(lambda) * (C / (eye(size(A)) - A/2/lambda));
Dd = 1/2/lambda * (C / (eye(size(A)) - A/2/lambda) * B) + D;

sysd = ss(Ad,Bd,Cd,Dd,Ts);

end

[John]

Thanks @Paul. That seems a useful approach.

I'm having some trouble with the code you posted, though. Although the step response looks the same, the two SS are different, and Simulink doesn't work with the manual version (unstable or other odd responses like very low gain).

Direct c2d(plant, Ts, 'tustin'):

A = 
        x1         x2
x1     0.9999  1.994e-05
x2     -5.037     0.9936

B = 
        u
x1  5.037e-05
x2      5.037

C = 
                x       x2
theta          1  9.968e-06

D = 
u
theta  2.518e-05

Sample time: 2e-05 seconds
Discrete-time state-space model.

And using the manual bilinear:

A = 
              x1         x2
   x1     0.9999  1.994e-05
   x2     -5.037     0.9936
 
  B = 
            u1
   x1  0.01126
   x2     1126
 
  C = 
              x1         x2
   y1   0.004472  4.458e-08
 
  D = 
              u1
   y1  2.518e-05
 
Sample time: 2e-05 seconds
Discrete-time state-space model.

So I think what's missing at the end of your function is a normalization of B, C, so that C(1) = 1.

So C --> C * C(1)

So B --> B / C(1)

Otherwise the responses (if broken out into matrices) will not match. Any thoughts on this?

[Paul]

I can't comment on your example because I don't know its original continuous time plant model. However, I poked around in c2d and the only diference between it and tempbilinear is the scaling of the Bd and Cd matrices.

Here's an example

rng(100);
sysc = rss(3,1,1);
Ts = 0.05;
sysd1 = c2d(sysc,Ts,'tustin');
sysd2 = tempbilinear(sysc,Ts);

Verify that sysd1 and sysd2 have the same input/output response.

zpk(sysd1)

ans =
 
  0.028499 (z+1) (z^2 - 1.806z + 0.8158)
  --------------------------------------
     (z-0.8195) (z-0.9464) (z-0.9802)
 
Sample time: 0.05 seconds
Discrete-time zero/pole/gain model.

zpk(sysd2)

ans =
 
  0.028499 (z+1) (z^2 - 1.806z + 0.8158)
  --------------------------------------
     (z-0.8195) (z-0.9464) (z-0.9802)
 
Sample time: 0.05 seconds
Discrete-time zero/pole/gain model.

Comparing the state space model elements

sysd2.a./sysd1.a  % A is the same

ans = 3×3

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

sysd2.b./sysd1.b

ans = 3×1

4.4721 4.4721 4.4721

sysd2.c./sysd1.c

ans = 1×3

0.2236 0.2236 0.2236

sysd2.d./sysd1.d  % D is the same

ans = 1.0000

The c2d result can be recovered by applying scale factors of

[sqrt(Ts) 1./sqrt(Ts)]

ans = 1×2

0.2236 4.4721

to the b and c matrices of sysd2 (or just change the code in tempbilinear), if it's important that tempbilinear yield the same state trajectories as c2d.

Alternatively, include the state variables in the output vector of sys, and then use the corresponding outputs of sysd2 so as to not be concerned with the internal dynamics of sysd2.

I don't know what's going on in Simulink, but everything should be fine using sysd1 or sysd2 as long as the only their output is used (because they have the same transfer function). If using the state vector for feedback, then, yes, there will be differences, and those can be addressed by a simple modification to tempbilnear if the c2d state response is preferred, or inluding the states in the output vector.

Also, I did some experimenting and found that applying the uncertainty to the continuous uss and then applying c2d is equivalent (in an I/O sense) to applying tempbilinear to the continuous uss and then applying the same uncertainty to the discrete uss.

Ts = 1/50e3;
% Any values are a stable system with expected 2nd-O response
%P = {};
P.J_kgPm2 = 1.0000e-05;
P.b = 3.2000e-03;
P.k_kgPm = 2.5266e+00;
multJ_u = ureal('multJ_u', 1, 'percent', 10);
multB_u = ureal('multB_u', 1, 'percent', 5);
multK_u = ureal('multK_u', 1, 'percent', 5);
%P_u = {};
P_u.J_kgPm2 = P.J_kgPm2 * multJ_u;
P_u.b = P.b * multB_u;
P_u.k_kgPm = P.k_kgPm * multK_u;
% Make continuous-time uss system
A = [0, 1; -P_u.k_kgPm / P_u.J_kgPm2, -P_u.b / P_u.J_kgPm2];
B = [0; 1 / P_u.J_kgPm2];
C = [1 0];
D = [0];
sys = ss(A, B, C, D);
sys = ss(A, B/dcgain(sys), C, D, ...
    'StateName', {'theta', 'angVel'}, ...
    'InputName', {'u'}, ...
    'OutputName', {'theta'});

Generate 100 samples of sys and discretize them

numModels = 100; % reduced to meet Answers runtime requirement
[csys,samplevalues] = usample(sys,numModels);
sysd1_u = c2d(csys, Ts, 'tustin');

Generate the discretized uncertain model

sysd = tempbilinear(sys,Ts);

Apply the uncertainty samples to the discrete-time uss model

sysd2_u = usubs(sysd,samplevalues);

Compare. The frequency response of their difference is essentially zero.

bode(sysd1_u-sysd2_u)

function sysd = tempbilinear(sysc,Ts)
A = sysc.A;
B = sysc.B;
C = sysc.C;
D = sysc.D;
lambda = 1/Ts;
Ad = (eye(size(A)) - A/2/lambda) \ (eye(size(A)) + A/2/lambda);
Bd = 1/sqrt(lambda) * ((eye(size(A)) - A/2/lambda) \ B);
Cd = 1/sqrt(lambda) * (C / (eye(size(A)) - A/2/lambda));
Dd = 1/2/lambda * (C / (eye(size(A)) - A/2/lambda) * B) + D;
sysd = ss(Ad,Bd,Cd,Dd,Ts);
end