Step() not working as expected with USS state space (Uncertain Sys)
8 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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')
0 comentarios
Respuesta aceptada
Paul
el 26 de Abr. de 2023
Editada: Paul
el 3 de Mayo de 2023
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);
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);
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 comentarios
Paul
el 28 de Abr. de 2023
Editada: Paul
el 30 de Abr. de 2023
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)
zpk(sysd2)
Comparing the state space model elements
sysd2.a./sysd1.a % A is the same
sysd2.b./sysd1.b
sysd2.c./sysd1.c
sysd2.d./sysd1.d % D is the same
The c2d result can be recovered by applying scale factors of
[sqrt(Ts) 1./sqrt(Ts)]
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
Más respuestas (0)
Ver también
Categorías
Más información sobre Linear Algebra en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!