FRD MIMO plant closed loop control, incorrect response

2 views (last 30 days)
  • I have an identified MIMO plant with 3 inputs and 4 outputs. I use a MIMO FRD object to store the frequency response of the system.
  • I am creating a control loop on 2 of the outputs of the system. I use the matlab function feedback to close the loop between the controller (it is an PPF controller) and the plant. See the first figure for the controller, and the second figure for the plant, open loop and closed loop responses.
  • The controlled input/output plots (plots [3,2] and [4,3]) seem to give a correct response, as well as the directly connected parts of the MIMO system (all plots except [1,1] and [2,1]) but the controller seems to have no effect on the indirectly connected parts ([z_inp->z1] and [z_inp -> z2]).
  • This seems to be because it is a numerical (frd) model. If I generate a similar system with analytical transfer functions, the controller does have an effect. The following figure shows a (badly) analytically identified system with tfest which is controlled using the exact same controller. notice that the plots for [z_inp->z1] and [z_inp->z2] are effected by the controllers that control the plots [pzt_u1->pzt_y1] and [pzt_u2->pzt_y2], which is expected.
Why is the frd model not effected correctly by the controllers, and can I fix this so the numerical model shows the correct response?
Below is the relevant code. If you need more info please ask, but I cannot share the entire code because of intellectual property reasons.
% some constants
fs = 5000;
ts = 1/fs;
maxfreq = 2e2;
f1 = logspace(0,log10(maxfreq),fs);
window_1 = barthannwin(2^15);
win_overlap = []; %default is [] , which is 50% overlap
% I do the following for each transfer function. the only change is the
% input/output names and the data (u31,y31)
T31 = tfestimate(u31, y31, window_1,win_overlap,f1,fs);
T31 = frd(T31,2*pi*f1,ts, 'u', 'z_inp', 'y', 'pzt_y1');
%consequently I combine the frd models into a MIMO model:
tfs = fcat([ T11, T12, T13; ...
T21, T22, T23; ...
T31, T32, T33; ...
T41, T42, T43 ]);
% make controller
g = 0.5;
wp = 76.9*2*pi;
d = 0.3;
den = [1/wp^2 2*d/wp 1];
Ctot = tf({0 0 0 0; 0 0 -g 0; 0 0 0 -g},{1 1 1 1; 1 1 den 1; 1 1 1 den});
Ctot.y = {'z_inp', 'pzt_u1', 'pzt_u2'};
Ctot.u = {'z1', 'z2', 'pzt_y1', 'pzt_y2'};
Ctotd = c2d(Ctot, ts); %discretize controller
% close the loop
Ttot = feedback(tfs, Ctotd, -1, 'name');
Note: I use standard plot notation [y,x] where y is the row index and x is the column index, starting from the top left.
  3 Comments
Castor Verhoog
Castor Verhoog on 24 Jan 2023
Edited: Castor Verhoog on 24 Jan 2023
@Paul Thank you for your comment. Interesting results, I would indeed expect the frd system to behave the same or similar as the ss system. Of course, a main difference between your 1 input 2 output system is that it is not actually MIMO (you could call it SIMO). I attached a .mat file to my original question with all frd models Tij. You can merge them yourself and see if you can duplicate my results.
Paul
Paul on 24 Jan 2023
I'll take a look when I get chance. Not sure what you mean by "1 input 2 output system." Ttot has 4 outputs and 3 inputs. I was only showing the response of (1,1) and (2,1) elements becuase those were the ones of concern in the Question.

Sign in to comment.

Answers (1)

Paul
Paul on 25 Jan 2023
Edited: Paul on 25 Jan 2023
Define the controller
fs = 5000;
ts = 1/fs;
g = 0.5;
wp = 76.9*2*pi;
d = 0.3;
den = [1/wp^2 2*d/wp 1];
Ctot = tf({0 0 0 0; 0 0 -g 0; 0 0 0 -g},{1 1 1 1; 1 1 den 1; 1 1 1 den});
Ctot.y = {'z_inp', 'pzt_u1', 'pzt_u2'};
Ctot.u = {'z1', 'z2', 'pzt_y1', 'pzt_y2'};
Ctotd = c2d(Ctot, ts); % discretize controller
Load the plant transfer functions
load(websave('tfunctions','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1272540/tfunctions.mat'))
Verify that TFs have already been converted to discrete time FRD LTI objects with sampling period = 1/5000 and the final point in the frequency vector corresponds to 200 Hz
T11.FrequencyUnit
ans = 'rad/TimeUnit'
T11.Ts
ans = 2.0000e-04
T11.Frequency(end)/2/pi
ans = 200.0000
Create the plant, use P instead of tfs
P = [ T11, T12, T13; ...
T21, T22, T23; ...
T31, T32, T33; ...
T41, T42, T43 ];
P.InputName = Ctot.y;
P.OutputName = Ctot.u;
Plot P(1,1), it looks like the corresponding plot in the original question
opts = bodeoptions;
opts.FreqUnits = 'Hz';
opts.XLim = {[10 200]};
bodeplot(T11,opts)
Form closed-loop system
Ttot = feedback(P, Ctotd, 'name', -1);
Compare the (1,1) terms of the plant and the closed-loop system
bodeplot(P(1,1),Ttot(1,1),opts)
They look the same, but let's see how close they really are
bodeplot(P(1,1)-Ttot(1,1),opts)
So they're close, but not really identical. They differ by about 1e-5 - 1e-6 which is small, but not really that small when testing for equality.
We can derive a closed form expression for Ttot(1,1) and compute it.
c = Ctotd(2,3);
sigma_4 = P(1,3)*c + (-P(1,2)*P(3,3) + P(1,3)*P(3,2))*c^2;
Warning: Ignoring all input names because of name conflicts.
sigma_3 = P(1,2)*c + (+P(1,2)*P(4,3) - P(1,3)*P(4,2))*c^2;
Warning: Ignoring all input names because of name conflicts.
sigma_6 = (P(3,2) + P(4,3))*c + (P(3,2)*P(4,3) - P(3,3)*P(4,2))*c^2 + 1;
Warning: Ignoring all input names because of name conflicts.
Warning: Ignoring all output names because of name conflicts.
Warning: Ignoring all input names because of name conflicts.
D11 = (P(4,1)*sigma_4 + P(3,1)*sigma_3)/sigma_6;
Warning: Ignoring all output names because of name conflicts.
Ttot11 = P(1,1) - D11;
Warning: Ignoring all input names because of name conflicts.
Verify that Ttot11 is numerically equivalent to Ttot(1,1). The difference is about 1e-15 - 1e-16, so to within numerical accuracy.
bodeplot(Ttot11 - Ttot(1,1),opts)
From the form of Ttot11, it must be the case that D11 is the difference between P(1,1) and Ttot(1,1)
bodeplot(D11,opts)
This plot is identical (or at least nearly so) to the Bode plot of the difference shown above.
In summary, the feedback does have a small effect on the (1,1) term of the closed loop system. The effect is small because of how the plant TFs and the controller combine into D11. If needed, you can explore further and try to determine why D11 is that small and if that makes sense for your application.
I believe the same effect will apply for Ttot(2,1) as well.
Of course, all of the discussion above was limited to the frequency range 10-200 Hz. If we look to lower frequencies, the difference between P(1,1) and Ttot(1,1) is noticeable
opts.XLim = {[1 200]};
bodeplot(P(1,1),Ttot(1,1),opts)
Here is the derivation of the closed-loop transfer function matrix
% plant
P = sym('P',[4 3])
P = 
% controller
C = zeros(3,4,'sym');
syms c
C(2,3) = c;
C(3,4) = c
C = 
% closed-loop
H = inv(eye(4) + P*C)*P
H = 

Community Treasure Hunt

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

Start Hunting!