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

2 次查看（过去 30 天）
John
John 2023-4-26
编辑： Paul 2023-5-3
Ran in:
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')

请先登录，再回答此问题。

采纳的回答

Paul
Paul 2023-4-26
编辑：Paul 2023-5-3
Ran in:
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;
  4 个评论
显示 2更早的评论
John
John 2023-4-27
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
Paul 2023-4-28
编辑：Paul 2023-4-30
Ran in:
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

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Data Extraction 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by