Convert discrete sys to continuous sys using delays using linear-fractional transformation (LFT)
3 次查看(过去 30 天)
显示 更早的评论
Hi,
I am trying to create a function that reutrns the continuous time MIMO transfer function of a discrete time sys. Howevee, for some reason although the bodeplot of the continuous system provides the same magnitude magnitude as the discrete sys, the phase angle seems to be mirrored at 0 deg. The function that I would like support is called d2c_exact and is given below.
Could someone figure out what I might be doing wrong? Thanks.
NOTE: I would like a solution where matrices operations are visible and not through matlab functions that obfuscate the process.
% EXAMPLE CODE %
% Declare MIMO matrices A,B,C,D
A = [-0.31261936441776167,0.28772503446055314;0.04760617923328707,0.8887042134058204];
B = [-0.607364659998142,-1.56763974476322,0.5900308503812629;1.4606427684802321,-1.6261530164028424,-0.8222697125090287];
C = [-0.004231925948322231,0.6979760136195045;-0.6958350856345005,-1.5121349753857836];
D = [-1.9523633596883074,-1.4728789255615329,0.14037581452786482;0.22729522014518927,1.2609066299547338,2.830197741550703];
% Define delay of 1 s
Ts = 1;
sys = ss(A,B,C,D,Ts);
% Discrete system
Pd = sys;
% Function that is the topic of the question
Pdc = d2c_exact(sys);
% Matlab equivalent function
Pc = d2c(sys);
% Create plots
wd = logspace(-2,2,200);
opts = bodeoptions;
opts.FreqUnits = 'rad/s';
opts.FreqScale = 'log';
opts.MagUnits = 'abs';
opts.MagScale = 'log';
opts.PhaseUnits = 'deg';
opts.Grid = 'on';
figure
bodeplot(Pd,wd,'b',opts)
hold on
bodeplot(Pdc,wd,'y',opts)
bodeplot(Pc,wd,'r--',opts)
hold off
h = get(gcf,'Children');
xline(h(2),fs/2,'k--');
xline(h(3),fs/2,'k--');
legend('Pd(z)','Pd(s)*inv(ZOH(s))','Pd(z) (ZOH sampling)','Nyquist freq. (fs/2)')
function sysc = d2c_exact(sysd)
nx = size(sysd.A,2);
s=tf('s');
Ts = sysd.Ts;
[A,B,C,D] = ssdata(sysd);
ZOH = exp(-Ts*s); % Transfer function of delay
M = append(ZOH,ZOH);
LR = M - ss(A + eye(nx));
sysc = C*feedback(eye(nx), LR)*B + D;
end
2 个评论
Paul
2023-12-20
Hi Joan,
The frqeuency response of Pdc is the conjugate of the others.
Do you have a source for the code in d2c_exact?
采纳的回答
Paul
2023-12-20
编辑:Paul
2023-12-20
Hi Joan,
I reviewed the link from the julia site and I think I know what they're doing. But I'm not 100% sure because it seems like they've made the SISO tf solution in post #8 more complicated than it needs to be, and perhaps also in post #9 (which is your link).
As best I can tell, the julia solution can't be implemented in Matlab (I'll show why below), but I think we can do something equivalent.
Here's the code with plots and my comments.
% EXAMPLE CODE %
% Declare MIMO matrices A,B,C,D
A = [-0.31261936441776167,0.28772503446055314;0.04760617923328707,0.8887042134058204];
B = [-0.607364659998142,-1.56763974476322,0.5900308503812629;1.4606427684802321,-1.6261530164028424,-0.8222697125090287];
C = [-0.004231925948322231,0.6979760136195045;-0.6958350856345005,-1.5121349753857836];
D = [-1.9523633596883074,-1.4728789255615329,0.14037581452786482;0.22729522014518927,1.2609066299547338,2.830197741550703];
% Define delay of 1 s
Ts = 1;
sys = ss(A,B,C,D,Ts);
% Discrete system
Pd = sys;
% Function that is the topic of the question
Pdc = d2c_exact(sys);
This isn't really the equivalent function. d2c by default assumes that sys includes the affect of a ZOH on the input, but I don't think that's necessarily an assumption to make here.
% Matlab equivalent function
Pczoh = d2c(sys);
For this example, foh works much better
Pcfoh = d2c(sys,'foh');
These warning suggests that sys was not developed as discrete-time approxiation using some form of z = exp(s*T).
Here's the new function.
Pdcnew = d2cnew(sys);
Create plots, but only up to the Nyquist frequency, because that's all that makes sense for Pd. Good match, particularly between Pd and Pdcnew.
wd = logspace(-2,2,200);
opts = bodeoptions;
opts.PhaseWrapping = 'on';
bodeplot(Pd,Pdcnew,Pcfoh,wd(wd<pi/Ts),opts)
The model Pdcnew is a system with no states and only internal delays.
Pdcnew
Here is your implementation of the post #9 in that julia thread.
function sysc = d2c_exact(sysd)
nx = size(sysd.A,2);
s=tf('s');
Ts = sysd.Ts;
[A,B,C,D] = ssdata(sysd);
The comment here is correct, i.e., that this is the transfer function of a delay. But the variable name is misleading because this is not the transfer function of a zero-order-hold.
ZOH = exp(-Ts*s); % Transfer function of delay
But the real problem is that the julia code for that line is
z = delay(-T)
which means that z represents an advance of T seconds (not a delay), which is consistent with the meaning z in a z-transform sense. That means that the equivalent Matlab line would be
z = exp(s*Ts)
But that will cause an error in Matlab because the Control System Toolbox only models delays, not advances.
M = append(ZOH,ZOH);
LR = M - ss(A + eye(nx));
sysc = C*feedback(eye(nx), LR)*B + D;
end
As best I can tell, the julia code approach is to form the block diagram of the discrete-time state space model and then replace z with exp(sT), or equivalently replace z^-1 with exp(-sT), and then form the continuous-time model from the block diagram. I don't understand why they implemented it the way they did, perhaps there's some subtle point that I don't appreciate.
Anyway, here's how that approach (as best I understand it) works in the Control System Toolbox.
function sysc = d2cnew(sysd)
Ts = sysd.Ts;
[A,B,C,D] = ssdata(sysd);
invz = ss(eye(size(A)),'InputDelay',Ts);
%sysc = parallel(series(series(B,feedback(invz,A,+1)),C),D);
sysc = C*feedback(invz,A,+1)*B + D;
end
更多回答(1 个)
Sulaymon Eshkabilov
2023-12-20
Is this what you are getting:
% EXAMPLE CODE %
% Declare MIMO matrices A,B,C,D
A = [-0.31261936441776167,0.28772503446055314;0.04760617923328707,0.8887042134058204];
B = [-0.607364659998142,-1.56763974476322,0.5900308503812629;1.4606427684802321,-1.6261530164028424,-0.8222697125090287];
C = [-0.004231925948322231,0.6979760136195045;-0.6958350856345005,-1.5121349753857836];
D = [-1.9523633596883074,-1.4728789255615329,0.14037581452786482;0.22729522014518927,1.2609066299547338,2.830197741550703];
% Define delay of 1 s
Ts = 1;
sys = ss(A,B,C,D,Ts);
% Discrete system
Pd = sys;
% Function that is the topic of the question
Pdc = d2c_exact(sys);
% Matlab equivalent function
Pc = d2c(sys);
% Create plots
wd = logspace(-2,2,200);
% Sampling freq:
fs=1e3;
figure
OPTs = bodeoptions;
OPTs.Grid = 'on';
OPTs.FreqScale = 'linear';
OPTs.Title.String = 'Bode Plot of Transfer Function';
bodeplot(Pd,wd,'b',OPTs)
hold on
bodeplot(Pdc,wd,'y',OPTs)
bodeplot(Pc,wd,'r--',OPTs)
hold off
h = get(gcf,'Children');
xline(h(2),fs/2,'k--');
xline(h(3),fs/2,'k--');
legend('Pd(z)','Pd(s)*inv(ZOH(s))','Pd(z) (ZOH sampling)','Nyquist freq. (fs/2)')
function sysc = d2c_exact(sysd)
nx = size(sysd.A,2);
s=tf('s');
Ts = sysd.Ts;
[A,B,C,D] = ssdata(sysd);
ZOH = exp(-Ts*s); % Transfer function of delay
M = append(ZOH,ZOH);
LR = M - ss(A + eye(nx));
sysc = C*feedback(eye(nx), LR)*B + D;
end
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Uncertainty Analysis 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!