LQI with time-delay system

Matlab's documentation mentions that state-space systems with dead-time (eg. plant or measurement delays) can be defined with ss(), using 'inputDelay' and 'outputDelay'.

But, documentation (and Matlab error messages) also say that time-delay Pade approximations must be used for lqr() or lqi(), ie those functions cannot take ss() with time delays as defined above with 'inputDelay' or 'outputDelay'.

If I have a feedback loop period (the delay in the measurement I'm looking to take into account), I can use lqi() with a discretized model input, which should yield the correct gains for the corresponding loop delay.

But for a continuous system modeled in Simulink with a transport delay, there is no discretization step, and I'd need to find lqi() gains using the continuous system as approximated by the model and pade.

But in this case, what is the Q used? I'm unclear on what Matlab would like as Q matrix input to lqi() for a continuous system.

As a guess, I tried the following, which yielded error messages. Say the loop period is ts, and I have some basic 2nd-order plant model (attached).

load("lqiWithDelay")
plant
ts = 1 / 100e3;
[num, den] = pade(ts, 3); % Nth-order
loopDelayPadeApproxTf = tf(num, den);
plantWithDelay = plant * loopDelayPadeApproxTf;

% Find non-delayed CL gains
wnTarget_rPs = 10 * 2 * pi;
H_lqi = [2 * 0.7 * wnTarget_rPs, 0.5, -wnTarget_rPs^2];
Q_lqi = H_lqi.'* H_lqi;
R_lqi = 1;
K_lqi = lqi(model, Q_lqi, R_lqi)

% Find delayed CL gains
neededDelayedQsize = size(plantWithDelay.A, 1) + 1; % +1 needed for LQI
neededQlqiSize = size(modelWithDelay.A, 1) + 1
bigDelayedQ = zeros(neededDelayedQsize);
bigDelayedQ(1:3, 1:3) = Q_lqi
K_lqi_delayed = lqi(modelWithDelay, bigDelayedQ, R_lqi)

This uses zeros(), and i also tried eye() -- but im only guessing, since i don't know what lqi() does under the hood for those additional pade states that lead to a bigger needed Q matrix.

Answer 1

[Paul]

Yes, if using a Pade apporoximant for the time delay, then the Q matrix will have to be expanded to include the additional state(s). In principle, you could buld the new Q matrix using the symmetric root locus technique, as I think is being done here, to drive the closed-loop poles to desirable locations. Maybe not penalizing that state(s) might be reasonable as well.

In addition, the Pade state(s) isn't available for feedback in the actual system, so they'd have to be reconstructed somehow, perhaps with a reduced order state estimator.

In the question you mentioned that the time delay is on the plant measurements, but the code has the time delay on the plant input. I'm going to assume it's on the input for what follows.

For this particular problem, the time delay at the input to the plant seems to be small enough that it is easily accomodated by the very large stability margins at the plant input that are guaranteed by LQ state feedback.

load("lqiWithDelay")
plant;
ts = 1 / 100e3;
%[num, den] = pade(ts, 3); % Nth-order
%loopDelayPadeApproxTf = tf(num, den);
%plantWithDelay = plant * loopDelayPadeApproxTf;  % delay is at the plant input

% Find non-delayed CL gains
wnTarget_rPs = 10 * 2 * pi;
H_lqi = [2 * 0.7 * wnTarget_rPs, 0.5, -wnTarget_rPs^2];
Q_lqi = H_lqi.'* H_lqi;
R_lqi = 1;

%K_lqi = lqi(model, Q_lqi, R_lqi) % what was model? not defined ... assume it's supposed to be plant
K_lqi = lqi(plant, Q_lqi, R_lqi);

sysp = ss(plant.A,plant.B,[eye(2);plant.C],[zeros(2,1);plant.D],'InputName','u','OutputName',{'x1' 'x2' 'y'});
int = tf(K_lqi(3),[1 0],'InputName','e','OutputName','ui');
Kfb = ss(K_lqi(1:2),'InputName',{'x1' 'x2'},'OutputName','ufb');
esum = sumblk('e = r - y');
usum = sumblk('u = -ui - ufb');

sysc = connect(sysp,int,Kfb,esum,usum,'r','y','u');

Bode plot of the open-loop transfer function (w/o delay) at the plant input and of the delay.

L = getLoopTransfer(sysc,'u',-1);
figure
bode(L,exp(-tf('s')*ts))

As expected from an LQ design, the gain margin is infinite and the phase margin is quite large, about 90 deg in this case. The phase lag that would be introduced by the time delay at the input to the plant is quite small at the open-loop gain crossover frequency (~5000 rads/sec) and the phase margin will hardly be changed. The additional phase lag will cause an open-loop phase cross-over at very high frequency, but there the gain will be quite low and the gain margin quite large.

The effect of the delay on the closed-loop step response is negligible.

sysp = ss(plant.A,plant.B,[eye(2);plant.C],[zeros(2,1);plant.D],'InputName','u','OutputName',{'x1' 'x2' 'y'},'InputDelay',ts);
int = tf(K_lqi(3),[1 0],'InputName','e','OutputName','ui');
Kfb = ss(K_lqi(1:2),'InputName',{'x1' 'x2'},'OutputName','ufb');
esum = sumblk('e = r - y');
usum = sumblk('u = -ui - ufb');

syscdelay = connect(sysp,int,Kfb,esum,usum,'r','y','u');

figure
step(sysc,syscdelay)

[John]

Thanks @Paul. Yes, you're partially correct, I was looking to do manually do a delay on the output, but as far as I know TF*exp() is the only method...since the input is delayed, as far as the measurement is concerned the output is also delayed -- there's no information about actual plant motion since there's no independent observer, so a delayed y results in either case. Is there a better way to represent measurement delay, aside from ss('outpuDelay'), which doesn'isn't compatible per my post?

" what was model? not defined ... assume it's supposed to be plant"

Thanks -- woops. Yes, as you surmised it was plant.

load("lqiWithDelay")
plant
ts = 1 / 100e3;
[num, den] = pade(ts, 3); % Nth-order
loopDelayPadeApproxTf = tf(num, den);
plantWithDelay = plant * loopDelayPadeApproxTf;
% Find non-delayed CL gains
wnTarget_rPs = 10 * 2 * pi;
H_lqi = [2 * 0.7 * wnTarget_rPs, 0.5, -wnTarget_rPs^2];
Q_lqi = H_lqi.'* H_lqi;
R_lqi = 1;
K_lqi = lqi(plant, Q_lqi, R_lqi)

In any case, agreed that should be inside the LQ margin, so let me pick another example (I'm just trying to get the methodology that Matlab will accept, vs an answer for whether the sys is stable or not). Check out the new attached lqiWithDelay, and my SL model with a 100 us delay. This is unstable, so i'd like to programmatically find the new gains via lqi() that would make it stable (without using phase margin).

From basic theory, it's clear that a delay T represents -w*T phase loss at that given frequency, so the phase margin is known once a CL system is designed (as you've also shown above). But, I'd like to know how to find gains to make an LQI loop stable with a given measurement delay, by finding the gains with a design that explicitly includes a plant with delay.

Since Matlab doesn't take ss('outputDelay) as an input to lqi(), i'm trying to find the method that's supported, probably something using pade() and an expanded Q...

[Paul]

InputDelay and OutputDelay can definitely have different implications depending on what one is doing.

Either type can be set explicitly during model construction

sys1 = ss(-1,1,1,0,'InputDelay',0.1)

sys1 =
 
  A = 
       x1
   x1  -1
 
  B = 
       u1
   x1   1
 
  C = 
       x1
   y1   1
 
  D = 
       u1
   y1   0
 
  Input delays (seconds): 0.1 
 
Continuous-time state-space model.

sys2 = ss(-1,1,1,0,'OutputDelay',0.2)

sys2 =
 
  A = 
       x1
   x1  -1
 
  B = 
       u1
   x1   1
 
  C = 
       x1
   y1   1
 
  D = 
       u1
   y1   0
 
  Output delays (seconds): 0.2 
 
Continuous-time state-space model.

Or, we can multiple by exp(-T*s) on the right for an InputDelay or the left for an OutputDelay

sys3 = ss(-1,1,1,0)*exp(-0.1*tf('s'))

sys3 =
 
  A = 
       x1
   x1  -1
 
  B = 
       u1
   x1   1
 
  C = 
       x1
   y1   1
 
  D = 
       u1
   y1   0
 
  Input delays (seconds): 0.1 
 
Continuous-time state-space model.

sys4 = exp(-0.1*tf('s'))*ss(-1,1,1,0)

sys4 =
 
  A = 
       x1
   x1  -1
 
  B = 
       u1
   x1   1
 
  C = 
       x1
   y1   1
 
  D = 
       u1
   y1   0
 
  Output delays (seconds): 0.1 
 
Continuous-time state-space model.

Your code multiplied by the Pade on the right, so I assumed it was an input delay.

I'm afraid I'm not aware of any particular methodoloogy for designing LQ controllers for systems with delays. Have you tried a literature search?

[John]

Thanks Paul.

"Or, we can multiple by exp(-T*s) on the right for an InputDelay or the left for an OutputDelay"

Wait...how does order matter? I'm startled that input or output is controllable in this manner:D Isn't a matrix multiplied by an exponential the same, whether left-or-right multuplied? (Each matrix entry is multiplied by a transfer function, and that's commutative).

Or is this Matlab-specific notation? If so, where is it documented so that i can read more?

"Your code multiplied by the Pade on the right, so I assumed it was an input delay."

Yeah, makes sense since this is an officially supported method. I had no idea... I was just trying to create a delay that would be the same as far as the measurement was concerned, as though the measurement was delayed.

"I'm afraid I'm not aware of any particular methodoloogy for designing LQ controllers for systems with delays. Have you tried a literature search?"

Yes, I tried some literature searches, but nothing "normal". I had assumed (wrongly) that there was an easy way to represent a state space system with delay in Matlab, and input that into the various LQ() functions.

Since it's not clear how a larger Q (when augmented by Pade entries) would be constructed, i guess that path is unavailable (or not clear).

Thanks for the help above, and this thread can be left as-is :) If i end up finding anything in the future, i'll update the post for anyone else looking to do the same (I saw at least one other person asking a few years back about something similar).

[Paul]

The Control System Toolbox is being very precise in how it models delays. Here's the doc page,

Models with Time Delays - MATLAB & Simulink (mathworks.com)

I suggest reading at least the first few Topics.

With lti objects, the * operation is shorthand for series. For SISO transfer functions, order doesn't matter as far as the I/O response. But it does matter for MIMO transfer functions, because MIMO transfer functions, in general, do not commute because they are matrices. Mutiplying with by scalar exp(-T*s) also makes a difference with respect to where the delay is placed (input or outupt).

If the delay is the same on all the inputs, then the same I/O response can be achieved by putting the same delay on all of the outputs. However, the state variables for a realization of such systems will not be the same.

The location of the delay(s) will also have different affects on the closed loop I/O response as well if the delay is on the feedback but the system output is not delayed.

In your case, the plant has one input and three feedback signals (two states and the output to the integral control). If you can figure out how to use a Pade approximation for your LQI design, you would only need one Pade at the input, but three if on the feedback (one for each feedback signal). In the latter case, you'd have three times as many Pade states to deal with than the former, so it's different.

[John]

"The location of the delay(s) will also have different affects on the closed loop I/O response as well if the delay is on the feedback but the system output is not delayed."

That makes sesne, although it shouldn't have an impact on the closed-loop behavior of the plant.

"If you can figure out how to use a Pade approximation for your LQI design, you would only need one Pade at the input, but three if on the feedback (one for each feedback signal). "

Interesting...I'll look into it.

1) "...the open-loop transfer function (w/o delay) at the plant input and of the delay."

How did you know that 'u' was the correct analysis point for this? Is it the case that any common point (eg plant output) would be also suitable? Ie it seems that should be okay, as long as all FB loops are captured by the analysis point, to show the impact of the FB gain system.

2) "As expected from an LQ design, the gain margin is infinite and the phase margin is quite large, about 90 deg in this case. The phase lag that would be introduced by the time delay at the input to the plant is quite small at the open-loop gain crossover frequency (~5000 rads/sec) and the phase margin will hardly be changed. The additional phase lag will cause an open-loop phase cross-over at very high frequency, but there the gain will be quite low and the gain margin quite large.

The effect of the delay on the closed-loop step response is negligible."

Hmmm... I made a quick SL check on this, and the loop goes unstable with a 100 us delay (1/10e3), which is odd considering that 100 us would only introduces -26 deg phase at 5000 rad/s: -w*delay = - 5000 * -100e-6 = -0.5 rad ~= -29 deg. Ie, phase goes from -90 to -119 deg. Why would that go unstable?

Run the SL model, with the gains found here:

load("lqiWithDelay")
% Find non-delayed CL gains
wnTarget_rPs = 10 * 2 * pi; % Low 10 Hz
H_lqi = [2 * 0.7 * wnTarget_rPs, 0.5, -wnTarget_rPs^2];
Q_lqi = H_lqi.'* H_lqi;
R_lqi = 1;
K_lqi = lqi(plant, Q_lqi, R_lqi)

[Paul]

What I meant by the "location of the delay(s) will also have different affects" is that if G(s) is the plant and D(s) is exp(-s*T), then the closed loop transfer function with delay in in the forward path is

H(s) = G*(s)*D(s)/(1 + G(s)*D(s))

and if in the feeback path

H(s) = G*(s)/(1 + G(s)*D(s))

These are not the same.

1) I did the analysis at the input to the plant because, when I wrote that answer, I thought that's where the delay should be placed because you right-multiplied by the delay in this line in your original code:

plantWithDelay = plant * loopDelayPadeApproxTf;

The stability margins that the LQ design guarantees are only guaranteed at the input to the plant (even for multiple inputs, which is not the case here). But there are not any guaranteed margins at the output of the plant. So the analysis point does (or can) make a difference.

2) I'm afraid I can't check your Simulink model. What you can do is add an open-loop linearization point in your Simulink model at the input to the plant, linearize the model (either sing the app or linearize), and then generate the Bode plot of the open-loop system (don't forget to multiply by -1). It should match the Bode plot I showed above if we're using the same plant model and feedback gains.

Yes, a 100 micro-sec delay at the plant input, which is 10x larger than the ts in the orignal question, will have a phase lag of ~29 deg at 5000 rad/sec and shouldn't cause instability, and seems like it does not.

If the Bode plot from the Simulink model matches what I showed above, or below for that matter, but the step response is going unstable, then there might be an issue with your solver settings.

load("lqiWithDelay")
plant;

Use higher value of ts, 100 microsec

%ts = 1 / 100e3;
ts = 100e-6;

%[num, den] = pade(ts, 3); % Nth-order
%loopDelayPadeApproxTf = tf(num, den);
%plantWithDelay = plant * loopDelayPadeApproxTf;  % delay is at the plant input
% Find non-delayed CL gains
wnTarget_rPs = 10 * 2 * pi;
H_lqi = [2 * 0.7 * wnTarget_rPs, 0.5, -wnTarget_rPs^2];
Q_lqi = H_lqi.'* H_lqi;
R_lqi = 1;
%K_lqi = lqi(model, Q_lqi, R_lqi) % what was model? not defined ... assume it's supposed to be plant
K_lqi = lqi(plant, Q_lqi, R_lqi);
sysp = ss(plant.A,plant.B,[eye(2);plant.C],[zeros(2,1);plant.D],'InputName','u','OutputName',{'x1' 'x2' 'y'});
int = tf(K_lqi(3),[1 0],'InputName','e','OutputName','ui');
Kfb = ss(K_lqi(1:2),'InputName',{'x1' 'x2'},'OutputName','ufb');
esum = sumblk('e = r - y');
usum = sumblk('u = -ui - ufb');
sysc = connect(sysp,int,Kfb,esum,usum,'r','y','u');
sysp = ss(plant.A,plant.B,[eye(2);plant.C],[zeros(2,1);plant.D],'InputName','u','OutputName',{'x1' 'x2' 'y'},'InputDelay',ts);
sysp.InputDelay

ans = 1.0000e-04

int = tf(K_lqi(3),[1 0],'InputName','e','OutputName','ui');
Kfb = ss(K_lqi(1:2),'InputName',{'x1' 'x2'},'OutputName','ufb');
esum = sumblk('e = r - y');
usum = sumblk('u = -ui - ufb');
syscdelay = connect(sysp,int,Kfb,esum,usum,'r','y','u');

L = getLoopTransfer(sysc,'u',-1);
Ldelay = getLoopTransfer(syscdelay,'u',-1);
figure
bode(L,Ldelay,{1e1 1e4}),grid

figure
step(sysc,syscdelay)

If you have trouble getting the Simulink model to match my analysis here, I can probably build a Simulink model and post it (assuming I can replicate my own results ;) )

[John]

Dang it, I did it again (responding in the wrong box). Thanks for moving...

1) "The stability margins that the LQ design guarantees are only guaranteed at the input to the plant (even for multiple inputs, which is not the case here). But there are not any guaranteed margins at the output of the plant. So the analysis point does (or can) make a difference."

Hmmm...I'll have to think through this more...

2) "I'm afraid I can't check your Simulink model. What you can do is add an open-loop linearization point in your Simulink model at the input to the plant, linearize the model"

Just to check, can't check as in my model has an error, or for policy reasons?

Either way, I increased the solver resolution to 1e-8, with no change in behavior; this solver resolution should be past the point where it needs to be.

3) I added an open loop output linearization point to u, and an input perturbation linearization point to y; is this the correct process and same points? It seems to be the right process at least from this video, and the same u point you're using with getLoopTransfer().

The bode that i see (with Model Linearizer) is different:

[Paul]

2. Personal policy

3. No, that's not correct. Delete both of those linearization points. Then, add one linearization point on the u_FB line, and make it's type 'looptransfer'; the symbol is shown on this doc page in the table for 'Loop transfer function'.

I have no idea how that Derivative block is implemented for linearization. Might be o.k., I don't know. But the Derivative block is strongly not recommended for time domain simulation. Instead, modify the plant block so it outputs y AND x2, and use x2 as the feedback into K_lqi(2).

[John]

"2. Personal policy"

Ah, got it -- apologies then for continually attaching .slx...

"Delete both of those linearization points. Then, add one linearization point on the u_FB line"

I see. The result is different than what you're seeing, mainly because the phase is 90 instead of -90 (though the gain also differs). I think my junction signs are correct given that you're also using usum = sumblk('u = -ui - ufb').

I checked plant and K_lqi, which should match what's in lqiWithDelay.mat:

K_lqi = 9.1242e+02 2.2417e-01 -3.9478e+05

zpk(plant) =

From input "u" to output "theta":

63165

-----------------------

(s^2 + 160s + 6.317e04)

When I add the 100 us delays in the loop, the bode plot is the same -- ie no change in the phase. So, the bode check I'm using (in Model Linearizer) seems wrong in my setup, inverted phase aside.

Here's the step response with no delay (immediately below),

and with 100 us delay (below):

[Paul]

No aplogies needed. Keep in mind that someone else reading this thread might take a look at the Simulink model.

Something has changed on your end. I'm attaching to this comment the .mat file that you provided with the original question and then running the code from the original question (both of which are what I've been using all along).

load("lqiWithDelay")

The plant is different than what you're showing immediately above.

zpk(plant)

ans =
 
  From input "u" to output "theta":
        8882.6
  ------------------
  (s^2 + 60s + 8883)
 
Continuous-time zero/pole/gain model.

%ts = 1 / 100e3;
%[num, den] = pade(ts, 3); % Nth-order
%loopDelayPadeApproxTf = tf(num, den);
%plantWithDelay = plant * loopDelayPadeApproxTf;
% Find non-delayed CL gains
wnTarget_rPs = 10 * 2 * pi;
H_lqi = [2 * 0.7 * wnTarget_rPs, 0.5, -wnTarget_rPs^2];
Q_lqi = H_lqi.'* H_lqi;
R_lqi = 1;
format short e
%K_lqi = lqi(model, Q_lqi, R_lqi)
K_lqi = lqi(plant, Q_lqi, R_lqi)

K_lqi = 1×3

1.0e+00 * 8.7839e+01 5.1269e-01 -3.9478e+03

And the gains are different than what you're showing immediately above.

Did you set the Pade order (for linearization) value in the Transport Delay block parameters to a positive integer? If not, that block linearizes to a gain of 1.

[John]

"Did you set the Pade order (for linearization) value in the Transport Delay block parameters to a positive integer? If not, that block linearizes to a gain of 1"

Ah...thanks. I had no idea that needed to be changed, but makes sense. The phase delta is as expected now, after setting to the appropriate value.

"The plant is different than what you're showing immediately above"

Huh -- i'm baffled by this. I run

load("lqiWithDelay")

zpk(plant)

and get what i showed above. Clearly, there's been a change on my end somehow that propagated to the .mat file i'm using.

Either way, for the plant and gains that you have, the gain plot matches the plot you showed -- but my bode phase still show +90 deg, instead of -90 deg phase. Is the loop transfer function point set correctly, in the screenshot i shared?

[Paul]

Yes, it looks like the linearization point is now set correctly.

As previously stated: "What you can do is add an open-loop linearization point in your Simulink model at the input to the plant, linearize the model (either sing the app or linearize), and then generate the Bode plot of the open-loop system (don't forget to multiply by -1)." (emphasis added).

The linearize function (or the linearization app) just computes the transfer function around the loop at the linearization point, call it O(s) (for open-loop). The characteristic equation of the system in terms of O(s) is

Delta(s) = 1 - O(s)

However, the Nyquist criterion is based on a characteristic equation of the form

Delta(s) = 1 + L(s).

To get the correct L(s) to check stability margins we need L(s) = -1*O(s), hence the multiplication -1, and which is also why I used -1 as the third argument to getLoopTransfer in an upstream comment where I showed the Bode plot.

Suppose we have a simpie, system that has a plant G(s) and unity, negative feedback. The closed loop system is:

H(s) = G(s)/(1 + G(s))

Based on the denominator, we'd make a Bode plot of G(s) as the open loop system. But, if you implemented that model in Simulink and added the linearization point on the signal line that is the input to the plant, the output of linearize would be -G(s) because of the negative sign at the sum junction that forms the error signal. We'd have to multiply the output of linearize by -1 before doing any stability margin analysis.

I now see that there is a version of getLoopTransfer that is a command line interface to working with Simulink models, but I'm not familiar with it.

The doc pages for either version of getLoopTransfer have some examples that make it clear that the multiply by -1 is needed to get the correct loop transfer function for bode, margin, etc.

[John]

(I deleted my previous post, since I realized I was confused about some things Connect() is doing, in your post)

Thanks @Paul, your explanation makes sense for the -1 needed for getLoopTransfer(). I remember you saying the -1, but wasn't sure why -- this clarifies, thanks.

I realized I was confused about some things Connect() is doing, from your post setting up the loop transfer function

1) Why is ss(sysc) showing x3 as a state, ie unknown? How come that isn't needed to be defined in the buildup of sysc? I would have thought that it's 'y' at the least, given

sysp = ss(plant.A,plant.B,[eye(2);plant.C],[zeros(2,1);plant.D],'InputName','u','OutputName',{'x1' 'x2' 'y'})

ss(sysc)
A = 
               x1          x2          x3
   x1           0           1           0
   x2  -1.373e+13  -1.281e+07   1.234e+11
   x3  -2.621e+05           0           0
 
  B = 
               r
   x1          0
   x2          0
   x3  2.621e+05
 
  C = 
      x1  x2  x3
   y   1   0   0
 
  D = 
      r
   y  0

2) How did you know what system to build up here, with the augmented C and D? I'm not sure how these SS system are compatible...but mainly I'm not sure why C, D are made to be (x1, x2, y)

sysp = ss(plant.A,plant.B,[eye(2);plant.C],[zeros(2,1);plant.D],'InputName','u','OutputName',{'x1' 'x2' 'y'})

3) What's the purpose of having all 3 (r, y, and u) in building up sysc?

connect(sysp, fromErrToUint, Kfb, esum, usum, 'r', 'y', 'u')

What is the 'u' for? The documentation page says Inputs and Outputs, but what does >1 Output do? I would have thought it creates a sys for each output (r-->y, and also r-->u), but that's not the case.

Answered my own question: the u is the analysis point as you've sad, and connect() knows it's not an output (but is an analysis point), since a group of outputs would be added as {'y', 'u'}, not independently as 'y', 'u'. So the last arg group is the set of analysis points.

[Paul]

1) It looks like we're not using the same system. I'm using lqiWithDelay.mat posted with the original question

load("lqiWithDelay")

The original plant has state names defined

plant.StateName

ans = 2x1 cell array
    {'theta' }
    {'angVel'}

wnTarget_rPs = 10 * 2 * pi;
H_lqi = [2 * 0.7 * wnTarget_rPs, 0.5, -wnTarget_rPs^2];
Q_lqi = H_lqi.'* H_lqi;
R_lqi = 1;
%K_lqi = lqi(model, Q_lqi, R_lqi) % what was model? not defined ... assume it's supposed to be plant
K_lqi = lqi(plant, Q_lqi, R_lqi);

But when sysp is defined, it's not given any state names.

sysp = ss(plant.A,plant.B,[eye(2);plant.C],[zeros(2,1);plant.D],'InputName','u','OutputName',{'x1' 'x2' 'y'});
sysp.StateName

ans = 2x1 cell array
    {0x0 char}
    {0x0 char}

int = tf(K_lqi(3),[1 0],'InputName','e','OutputName','ui');
Kfb = ss(K_lqi(1:2),'InputName',{'x1' 'x2'},'OutputName','ufb');
esum = sumblk('e = r - y');
usum = sumblk('u = -ui - ufb');
sysc = connect(sysp,int,Kfb,esum,usum,'r','y','u');

sysp doesn't have any state names, and int can't have a state name (it's a tranfser function), so connect returns a model without state names

sysc.StateName

ans = 3x1 cell array
    {0x0 char}
    {0x0 char}
    {0x0 char}

If we give sysp state names

sysp.StateName = {'one','two'};

Then those state names are carried through to the output of connect()

sysc = connect(sysp,int,Kfb,esum,usum,'r','y','u');
sysc.StateName

ans = 3x1 cell array
    {'one'   }
    {'two'   }
    {0x0 char}

connect() only cares about the names of the inputs and outputs to each of the systems being connected, which are all defined. The state names of the ss systems are window dressing as far as connect() is concerned.

2) I created sysp by augmenting the original plant with two additional outputs being the states of sysp, which are the same as the states of plant. I guess I could have called those outputs theta and angVel as in the original plant. The output names are arbitrary as long as they pair up with inputs to other systems to be connected, which are in this case are Kfb and esum, and/or defined as outputs of the connected system, as is the case with 'y'.

3) Correct.

[John]

Thanks Paul. "It looks like we're not using the same system. I'm using lqiWithDelay.mat posted with the original question"

Sorry about that; you're right, and the specific system wasn't my focus... So good catch :D I was using another version just to illustrate...

"If we give sysp state names sysp.StateName = {'one','two'};

Then those state names are carried through to the output of connect()"

"connect() only cares about the names of the inputs and outputs to each of the systems being connected, which are all defined. The state names of the ss systems are window dressing as far as connect() is concerned."

Ah, I see. Those both make sense, thanks.

"I created sysp by augmenting the original plant with two additional outputs being the states of sysp, which are the same as the states of plant. I guess I could have called those outputs theta and angVel as in the original plant. The output names are arbitrary as long as they pair up with inputs to other systems to be connected, which are in this case are Kfb and esum, and/or defined as outputs of the connected system, as is the case with 'y'."

I'm not sure I totally follow, but think i partially get it: you're saying that for connect() to work, all inputs and outputs have to be defined in the original system (sysp) for connect() to know what to connect to? If so, I thought connect() can be used by connecting various inputs and outputs, so why is a single augmented plant necessary? Ie something like

ss(plant.A,plant.B,plant.C, plant.D,'InputName','u','OutputName', 'y')

And then the state x1, x2 would be another connect() to the feedback gains K(1:2)

In any case, I'm having trouble understanding why my simulink models are not matching the stability of the loop transfer function results. For example, here is a bode() for an LQI plant and no delay, and plant and LQI with 1 us delay. I compute what phase delay is allowed so that the open loop xfer func does not exceed (say) 60 deg phase margin at the crossover freq (in this case it's 1.7 us), and then build an open loop xfer func with that delay -- just as you did. The expectation is that the delayed TF will show 60 deg phase margin at the crossover, and that is exactly what i see (ie phase at -120 deg, ie 60 deg phase margin).

Then, I expect that if running a simulink model with the same delay, the loop should be stable (since there's still 60 deg PM). So I tried an even smaller delay (1us), and see instability:

Here is what's running:

Is this an improper use of the Transport Delay? I would have thought this acts like a system delay. What might be the issue with this block? The solver is running at 1e-8 s fixed-step, with no change at e-9 or e-10, so the solver seems fine.

[Paul]

"If so, I thought connect() can be used by connecting various inputs and outputs, so why is a single augmented plant necessary? Ie something like

ss(plant.A,plant.B,plant.C, plant.D,'InputName','u','OutputName', 'y')

And then the state x1, x2 would be another connect() to the feedback gains K(1:2)"

You said it yourself. connect is linking up the systems by linking the elements of OutputName of each system to same-named elements of InputName for all of the other systems. The elements of StateName for any system don't come into play. That's just how it works, see the doc page for more info. That's why I formed sysp so that the state variables {'x1', 'x2'} were also in the output vector to link up as inputs to the Kfb system, and an additional output {'y'} to link as input to esum to form the error signal. Because y and x1 are the same quantity, I suppose I could have made sysp with only two outputs {'x2' , 'y'} and then set Kfb.InputName = {'y' , 'x2'}.

As for the simulink model, what does the step response look like with the transport delays set to zero (if allowed) or removed?

[John]

"That's just how it works, see the doc page for more info. That's why I formed sysp so that the state variables {'x1', 'x2'} were also in the output vector to link up as inputs to the Kfb system, and an additional output {'y'} to link as input to esum to form the error signal."

I see... I read the doc page but thought it'd be possible to (for example) do the following: the plant has output x1, x2. To this plant, connect another section that takes x1,x2, and C,D, and forms y -- similar to what an actual computation would do (y = Cx + D). Similarly, the output of the plant is already x1, x2, so that should be linkable already. I'll think it over what you wrote a bit more, since I believe the root issue is that I'm not fully understanding what connect() is actually doing with the sysp, ie the ss() matrices.

"As for the simulink model, what does the step response look like with the transport delays set to zero (if allowed) or removed?"

If the two Transport Delays are commented through, the step response looks great: exactly as expected. (The block diagram then looks like what i attached in this post).

[Paul]

Need to get on the same page. Attach a .mat file in a response to this comment that contains the plant and the K_lqi that you're currently working with. Pick a new name for the file, i.e., do NOT call it 'lqiWithDelay.mat'.

[John]

I apologize -- I made a mistake. It's good that you asked :) One of the gains was incorrect due to a typo. So, this all seems consistent now for the continuous-time case; thanks!

A new question came out of that: how would connect() be used with a discrete-time system? The doc page didn't seem to show any examples, but might have missed it. It seems it'd be the same in principle, but practically it's not clear.

Assume a discrete-time system: it's the same continuous-time plant, but a discrete-time loop (with associated discrete gains) in the feedback; so x1, x2, y are updated at a given loop rate, and K_d gains are correct discrete-time gains.

Using your connect() and getLoopTransfer()method, it's straightforward to add the discrete-time gains, but it seems that the loop time would need to be inserted into each loop element somehow as inline elements to y, x1, x2, so that connect() creates the correct system.

a) To modify x1, x2 , I thought something like this would be needed:

Kfb = ss(K_d(1:2),'InputName',{'x1' 'x2'},'OutputName','ufb', loopTime_s);

or

Kfb = ss(K_d(1:2), loopTime_s, 'InputName',{'x1' 'x2'},'OutputName','ufb');

But neither seems accepted syntax, since ss() accepts loop time only in the format of ss(A, B, C, D, Ts) per the documentation (not with eg gain K). How would I use K_d(1:2) with discrete-time ss()?

b) To modify y, it's not clear where disc-time could be specified, since previously the y feedback wasn't explicitely called, as y was an output from sysp. The closest was:

int = tf(K_d(3),[1 0],'InputName','e','OutputName','ui');

How would connect() be used to make a discrete-time system with x1, x2, y being updated at a given loop rate?

[Paul]

Static gains don't have to be given a sample time when connecting with discrete-time systems.

Here's the continuous-time system.

load("lqiWithDelay")
% Find non-delayed CL gains
wnTarget_rPs = 10 * 2 * pi;
H_lqi = [2 * 0.7 * wnTarget_rPs, 0.5, -wnTarget_rPs^2];
Q_lqi = H_lqi.'* H_lqi;
R_lqi = 1;
K_lqi = lqi(plant, Q_lqi, R_lqi);
sysp = ss(plant.A,plant.B,[eye(2);plant.C],[zeros(2,1);plant.D],'InputName','u','OutputName',{'x1' 'x2' 'y'});
int = tf(K_lqi(3),[1 0],'InputName','e','OutputName','ui');
Kfb = ss(K_lqi(1:2),'InputName',{'x1' 'x2'},'OutputName','ufb');
esum = sumblk('e = r - y');
usum = sumblk('u = -ui - ufb');
sysc = connect(sysp,int,Kfb,esum,usum,'r','y','u');

Discretize the plant assuming a ZOH at the input

Ts = 0.0001;
syspdiscrete = c2d(sysp,Ts,'zoh');

c2d preserves the InputName and the OutputName

syspdiscrete.InputName

ans = 1x1 cell array
    {'u'}

Define a discrete-time integrator, multiplied by K_lqi(3)

intdiscrete = tf(Ts*K_lqi(3),[1 -1],Ts,'InputName','e','OutputName','ui');

Connect in the usual way and compare

syscdiscrete = connect(syspdiscrete,intdiscrete,Kfb,esum,usum,'r','y','u');

figure
step(sysc,syscdiscrete),legend('Cont','Disc')

You can add a sample time to a static gain ss (or any other) model by explicitly using the Ts property

Kfb = ss(K_lqi(1:2),'InputName',{'x1' 'x2'},'OutputName','ufb','Ts',Ts);

But it won't matter as far as connect() is concerned

syscdiscrete = connect(syspdiscrete,intdiscrete,Kfb,esum,usum,'r','y','u');

figure
step(sysc,syscdiscrete),legend('Cont','Disc')

When using connect() for a hybrid system modeled as a discrete-time system, it's up to you to first properly discretize all of the continuous subsystems before calling connect, which is also true when using other functions like feedback.

[John]

1) "Static gains don't have to be given a sample time when connecting with discrete-time systems."

Agreed; but I'm trying to show that this loop component has the rate transition in it, ie discrete-time operation, as opposed to the plant, which is cts-time.

"Discretize the plant assuming a ZOH at the input"

Would this be a different loop configuration then what i have, then? Ie the plant would operate in discrete time, which seems like it'd be a different system.

From my previous post:

"Assume a discrete-time system: it's the same continuous-time plant, but a discrete-time loop (with associated discrete gains) in the feedback; so x1, x2, y are updated at a given loop rate, and K_d gains are correct discrete-time gains."

Here's the block diagram, for clarity. It's the same cts-time plant, but discrete-time inputs are going to it (ie plant evolves in continuous time, in response to ZOH inputs). I'm trying to represent each branch as discrete-time, hence the Ts in K*u, and trying to add a Ts to the y feedback path.

Or would the the disc-time connect() example you showed represent the above? (Ignoring high-level behavior, which agreed will largely be equivalent. I'm looking to represent this system exactly as the block diagram shows it, though).

"You can add a sample time to a static gain ss (or any other) model by explicitly using the Ts property"

Ah, great, thanks. I didn't realize (might have missed it in the docs) that 'Ts' needed to be added.

2) I thought about your augmented plant explanation more, so I'll try to say it back to see if it's right:

A standard LTI state-space system represents x. = Ax + Bu, y = Cx+D. With your definitions for C [eye(2); (2 0], the output y = Cx + D seems like it'd not be 'y', but the vector [x1; x2; y] which includes 'y' itself. So in that sense, this isn't a state-space system that would be used for standard analysis/feedback, since y isn't necessarily sensible in the state-space sense; it's purely for Matlab's connect() requirements since x1, x2, y are all needed. But i wouldn't use sysp for anything real.

Is this close to what you were saying?

And if so, then in the specific FB loop setup of this example where x1 (theta), x2 (angvel), y(theta) are used, is sysp equivalent doing the following manualSysp? (changes are bolded)

justPlant = ss(plant.A, plant.B, plant.C, plant.D, 'InputName','u','OutputName','y')

Kfb = ss([K_lqi(1), K_lqi(2)*s],'InputName','y','OutputName','ufb');

and connect() them using

manualSysp = connect(plantNotAug, Kfb, 'InputName','u','OutputName',{'x1', 'x2', 'y')

[Paul]

1) AFAIK, the Control System Toolbox does not support hybrid systems. So functions like connect and feedback require that all of the input systems have the same sample time. Simulink supports time-domain simulation of hybrid systems. So, yes, the system that I produced from connect is not exactly same as what Simulink is simulating in the time domain. The former is an approximation of the latter. However, keep in mind that, when using Simulink's linearization functions for hybrid systems, those functions first convert all of the continuous time blocks to their discrete-time approximations (it does this very carefully, it does not convert each continuous-time block one at a time I don't believe) and then returns a discrete-time system as the result. You should see this in your model when you linearize to generate the Bode plot of the loop transfer function at the breakpoint on the u_fb signal. If the Simulink model contains more than one discrete sample-time, I believe the linearization functions find the slowest sample time, which might be a multiple of the sample times in the model, and will convert all of the discrete- and continuous-time blocks in the model to that sample time before returning the result. I'm not 100% certain, but I'm pretty sure that Simulink functions like linearize ultimately use connect after all the underlying block conversions, should any be necessary.

2) Yes, I defined sysp so that its output vector would have the three signals I needed to close loops using connect. As I said above, I could have defined sysp to only have two signals in its output vector. I don't know what you mean by "y isn't necessarily sensible in the state-space sense;" but yes, you get the point, which is that connect works by pairing up names in the InputName and OutputName properties of the subsystems to be connected. Also, the names of the input signals and the names of the output signals for the system that is generated by connect all have to be in InputName and OutputName properties of the subsystems.

However, this code won't work.

load lqiWithDelay.mat
wnTarget_rPs = 10 * 2 * pi;
H_lqi = [2 * 0.7 * wnTarget_rPs, 0.5, -wnTarget_rPs^2];
Q_lqi = H_lqi.'* H_lqi;
R_lqi = 1;
K_lqi = lqi(plant, Q_lqi, R_lqi);

Here, justPlant is an ss object with one input 'u' and one output 'y'

justPlant = ss(plant.A, plant.B, plant.C, plant.D, 'InputName','u','OutputName','y');
justPlant.InputName

ans = 1×1 cell array
    {'u'}

justPlant.OutputName

ans = 1×1 cell array
    {'y'}

which could have just as easily been accomplished by making them InputName and OutputName properties of plant. plant already has 'u' as the InputName

plant.InputName

ans = 1×1 cell array
    {'u'}

The OutputName is

plant.OutputName

ans = 1×1 cell array
    {'theta'}

which can be changed if you want

plant.OutputName = 'y';

Continuing with the code ....

These lines generates Kfb as ss with two inputs and one output

s = tf('s'); % I assume this is what you meant
% because s is a tf, so will be [K_lqi(1), K_lqi(2)*s], which has to be converted to ss first 
Kfb = ss([K_lqi(1), K_lqi(2)*s]);
% and then the InputName and OutputName can be assigned
Kfb.InputName = 'y';Kfb.OutputName = 'ufb';

Because only a single name was provided for InputName, but there are two inputs, the signal name was automagically expanded to two unique names.

Kfb.InputName

ans = 2×1 cell array
    {'y(1)'}
    {'y(2)'}

Note also that because Kfb was formed with a pure differentiator, it's actually in descriptor form

Kfb.E

ans = 2×2

0 1 0 0

I don't think that's a problem, but frankly I'm not sure how connect works if one or more subsystems are in descriptor form but the resulting system from connect doesn't need to be.

This line has multiple problems. Nowhere in this thread, nor on the doc page connect, did you see the InputName and OutputName name/value pairs. So this line will throw an error (even after adding the closing brace on the output name cell array and changing plantNotAug to justPlant?)

try
%    manualSysp = connect(plantNotAug, Kfb, 'InputName','u','OutputName',{'x1', 'x2', 'y'})
    manualSysp = connect(justPlant, Kfb, 'InputName','u','OutputName',{'x1', 'x2', 'y'});
catch ME
    ME.message
end

ans = 'Input argument 3 is not a dynamic system or has some unspecified I/O names.'

Even if we use the correct syntax, there will still be a problen because 'x1' and 'x2' are not OutputNames of either justPlant or Kfb, so connect doesn' know what you're asking for as the output signals from manualSysp.

try
    manualSysp = connect(justPlant, Kfb,'u',{'x1', 'x2', 'y'});
catch ME
    ME.message
end

ans = 'The signal name "x1" is not listed as output by any block.'

If we fix that problem,

try
    manualSysp = connect(justPlant, Kfb,'u','y');
catch ME
    ME.message
end

Warning: The following block inputs are not used: y(1),y(2).

Warning: The following block outputs are not used: ufb.

we don't get an error, but it's still wrong because the InputNames to Kfb are 'y(1)' and 'y(2)' and those names are not OutputNames from justPlant, and 'ufb', which is the OutputName from Kfb, is not the InputName to justPlant, which is 'u'.

So manualSysp is something

manualSysp

manualSysp =
 
  A = 
          x1     x2
   x1      0      1
   x2  -8883    -60
 
  B = 
          u
   x1     0
   x2  8883
 
  C = 
      x1  x2
   y   1   0
 
  D = 
      u
   y  0
 
  E = 
       x1  x2
   x1   1   0
   x2   0   1
 
Continuous-time state-space model.

but it's nothing that you want.

I suggest carefully reading through the examples on the connect doc page (don't worry about Index Based Interconnection), if you haven't already, then come back and see if what I did with connect in any of my previous comments makes more sense. It helps very much to draw the block diagram and carefully label every signal.

[John]

Thanks Paul. I'll answer (2) later, but looking at (1) and the response via Connect():

Discretizing the whole system does seem to work for lqiWithDelay.mat, so thanks for that example.

It doesn't seem to work for other system that are stable when there's only a discrete loop but continuous plant, so for some reason the methods don't seem to be interchangeable. For example, i try with this set of plant and gains (attached), and this has a stable discrete step in Simulink, but not in the step response with the Matlab connect() example (where the discrete version blows up)

Would you think it's the same to discretize the full FB system, vs using discrete gains?

Aside from a cts plant, that's the main difference i can think of.

Connect() results:

Just the cts portion:

Simulink results with the same params and block diagram from earlier, which is also at y=0.98 @t = 8e-3, so i don't think i have a typo this time :

[Paul]

load newPlantAndGains
Ts

Ts = 1.0000e-05

sysp = ss(plant.A,plant.B,[eye(2);plant.C],[zeros(2,1);plant.D],'InputName','u','OutputName',{'x1' 'x2' 'y'});
int = tf(K_lqi(3),[1 0],'InputName','e','OutputName','ui');
Kfb = ss(K_lqi(1:2),'InputName',{'x1' 'x2'},'OutputName','ufb');
esum = sumblk('e = r - y');
usum = sumblk('u = -ui - ufb');
sysc = connect(sysp,int,Kfb,esum,usum,'r','y','u');
figure
step(sysc);

Discretize the plant and define the discrete integrator with the Ts in the .mat file

syspdiscrete = c2d(sysp,Ts,'zoh');
intdiscrete  = tf(Ts*K_lqi(3),[1 -1],Ts,'InputName','e','OutputName','ui');
syscdiscrete = connect(syspdiscrete,intdiscrete,Kfb,esum,usum,'r','y','u');
figure
step(syscdiscrete)

Check the poles of the closed loop system

format short e
pole(sysc)

ans = 3×1

1.0e+00 * -2.8081e+05 -7.0370e+03 -5.0507e+02

The time constant of the fastest pole is

-1/ans(1)

ans = 
   3.5612e-06

Typically we'd want the sample time to be 5-10 faster than the fastest pole. Here we'll make it just a bit faster

Ts = 1e-6;

Redo with this faster sample time

syspdiscrete = c2d(sysp,Ts,'zoh');
intdiscrete  = tf(Ts*K_lqi(3),[1 -1],Ts,'InputName','e','OutputName','ui');
syscdiscrete = connect(syspdiscrete,intdiscrete,Kfb,esum,usum,'r','y','u');
figure
step(sysc,syscdiscrete)

[John]

Good point; but why does this work in Simulink, then? It seems the two systems should behave ~identically.

Maybe this difference relates to the plant being discrete vs continuous...ie maybe need to use the mixed-signal capability of the open loop xfer function that's available in Simulink (vs using connect() and an all-discrete matlab system), since the poles of the cts-plant version would be different.

To your point, that'd need the -1x for the simulink OLTF given how it computes, so perhaps needing to use the matlab command interface for the simulink version.

Edit: It's also not clear to me why changing wnTarget_rPs to very small doesn't seem to cause the discrete system's open loop xfer function around 'u', to ever cross 0 dB (eg using the plant and same sampling rate as in my last mat file). Of course, due to loop sampling rate, there is a max bode() limit at rate/2, and i would thought a very low control BW would allow the bode magnitude to cross 0 before that rate/2 limit: the high-freq of the loop xfer function includethe phase and noise behavior. But, that 0-crossing doesn't seem to occur even for wnTarget_rPs ~= plant resonance(ie very low control BW).

Is getLoopTransfer() different for discrete system? As before, in the docs i don't see any reference to using this function on discrete systems, so perhaps this is the wrong usage. I thought i could just do:

getLoopTransfer(syscdiscrete,'u',-1);

[Paul]

The two systems aren't the same. In Simulink the plant is integrated forward using the continuous-time solver, whereas in Matlab with the connect the plant is discretized. They can be close, but they're not the same.

Regardless, I didn't see it work in Simulink. I tried plant and K_lqi from newPlantAndGains.mat (plant is the same as what we've been using all along) in Simulink with a fixed step solver with step size of 1e-8 and a discrete-time integrator and ZOH blocks on the feedback paths, all with sample time Ts = 1e-5, the response went unstable and the simulation stopped with an error. When I set Ts = 1e-6, all was fine. I can probably post code here to recreate that Simulink result.

As previously discussed, the way Simulink linearizes a hybrid system is by linearizing the nonlinear blocks and then discretizing the continuous-time blocks and then connecting it all together, which is exactly what I did with connect. I expect using linearize on the model will return the same result as I got with connect. Maybe I'll try this.

Like (all?) other Control System Toolbox functions, getLoopTransfer works just fine with discrete-time sytems.

Demonstration:

load newPlantAndGains

sysp = ss(plant.A,plant.B,[eye(2);plant.C],[zeros(2,1);plant.D],'InputName','u','OutputName',{'x1' 'x2' 'y'});
int = tf(K_lqi(3),[1 0],'InputName','e','OutputName','ui');
Kfb = ss(K_lqi(1:2),'InputName',{'x1' 'x2'},'OutputName','ufb');
esum = sumblk('e = r - y');
usum = sumblk('u = -ui - ufb');
sysc = connect(sysp,int,Kfb,esum,usum,'r','y','u');

Ts = 1e-6;
syspdiscrete = c2d(sysp,Ts,'zoh');
intdiscrete  = tf(Ts*K_lqi(3),[1 -1],Ts,'InputName','e','OutputName','ui');
syscdiscrete = connect(syspdiscrete,intdiscrete,Kfb,esum,usum,'r','y','u');

figure
bode(getLoopTransfer(sysc,'u',-1),getLoopTransfer(syscdiscrete,'u',-1),{10 1e7});
grid

[John]

Thanks @Paul. I think there are two issues at play that are causing this difference:

• Per the discrete-time block diagram i posted, I was using Transport Delay after the rate transition -- but this doesn't work (i had assumed it would account for discrete-time, but it doesn't). It needs to be before the rate transition, or i need to use z^-n delay. Using Transport Delay in this manner my system unstable. Fixing this made the SL simulations at least somewhat match at low loop rates.
• There seems to be some difference between connect() and SL that i've been unable to understand so far, perhaps related to the discretization method you pointed out:

"In Simulink the plant is integrated forward using the continuous-time solver, whereas in Matlab with the connect the plant is discretized. They can be close, but they're not the same. ... As previously discussed, the way Simulink linearizes a hybrid system is by linearizing the nonlinear blocks and then discretizing the continuous-time blocks and then connecting it all together, which is exactly what I did with connect. I expect using linearize on the model will return the same result as I got with connect. Maybe I'll try this."

I'll pause on this post and start another specific thread on discrete time sims, that I'm unable to make even closely match the connect() versions, as that seems to be the limiting question at the moment...

[Paul]

In the Simulink model in this comment it looked like the Transport Delay blocks were Commented Through. Yes, the transport delay should go before the Rate Transition blocks.

I don't see any difference at all between the results from connect/getLoopTransfer and from linearizing the Simulink model using linearize (which I'm sure is the same as what you'd get using any other method to linearize the model).

Load the system parameters and set the sample time.

load newPlantAndGains

plant.C = [eye(2);plant.C];
plant.Ts

ans = 0

Ts = 1e-6;

Buld the Simulink model with the continuous-time plant, discrete-time integrator, and feedback signals sampled at Ts.

bdclose('all');
hsys = new_system('lqitest');

hstep = add_block('simulink/Sources/Constant','lqitest/Step');

herror = add_block('simulink/Math Operations/Subtract','lqitest/error');

%hint = add_block('simulink/Continuous/Integrator','lqitest/int');
hint = add_block('simulink/Discrete/Discrete-Time Integrator','lqitest/int');
set_param(hint,'SampleTime','Ts');

hintgain = add_block('simulink/Math Operations/Gain','lqitest/intgain');
set_param(hintgain,'Gain','K_lqi(3)')

hcontrol = add_block('simulink/Math Operations/Subtract','lqitest/control');
set_param(hcontrol,'Inputs','--');

hplant = add_block('cstblocks/LTI System','lqitest/plant');
set_param(hplant,'MaskValues',{'plant';'[]';'0'});

hdemux = add_block('simulink/Signal Routing/Demux','lqitest/demux');
set_param(hdemux,'Outputs','[2 1]');

hstategain = add_block('simulink/Math Operations/Gain','lqitest/stategain');
set_param(hstategain,'Gain','K_lqi(1:2)');
set_param(hstategain,'Multiplication','Matrix(K*u)');

% hdelay1 = add_block('simulink/Continuous/Transport Delay','lqitest/delay1');
% set_param(hdelay1,'DelayTime','0');
% 
% hdelay2 = add_block('simulink/Continuous/Transport Delay','lqitest/delay2');
% set_param(hdelay2,'DelayTime','0');

hzoh1 = add_block('simulink/Discrete/Zero-Order Hold','lqitest/zoh1');
set_param(hzoh1,'SampleTime','Ts');

hzoh2 = add_block('simulink/Discrete/Zero-Order Hold','lqitest/zoh2');
set_param(hzoh2,'SampleTime','Ts');

hout = add_block('simulink/Sinks/To Workspace','lqitest/output');
set_param(hout,'VariableName','y');

PH = 'PortHandles';
add_line(hsys, ...
    get_param(hstep,PH).Outport,...
    get_param(herror,PH).Inport(1));

add_line(hsys,...
    get_param(herror,PH).Outport,...
    get_param(hint,PH).Inport);

add_line(hsys,...
    get_param(hint,PH).Outport,...
    get_param(hintgain,PH).Inport);

add_line(hsys,...
    get_param(hintgain,PH).Outport,...
    get_param(hcontrol,PH).Inport(2));

add_line(hsys,...
    get_param(hstategain,PH).Outport,...
    get_param(hcontrol,PH).Inport(1));

add_line(hsys,...
    get_param(hcontrol,PH).Outport,...
    get_param(hplant,PH).Inport);

add_line(hsys,...
    get_param(hplant,PH).Outport,...
    get_param(hdemux,PH).Inport);

add_line(hsys,...
    get_param(hdemux,PH).Outport(2),...
    get_param(hout,PH).Inport);

add_line(hsys,...
    get_param(hdemux,PH).Outport(2),...
    get_param(hzoh2,PH).Inport);

add_line(hsys,...
    get_param(hdemux,PH).Outport(1),...
    get_param(hzoh1,PH).Inport);

add_line(hsys,...
    get_param(hzoh1,PH).Outport,...
    get_param(hstategain,PH).Inport);

add_line(hsys,...
    get_param(hzoh2,PH).Outport,...
    get_param(herror,PH).Inport(2));

Uncomment this line to open the model on the Simulink canvas for by-hand editing, execution, etc.

% open_system(hsys);

Simulate the model in the time domain with a fixed step solver with a simulation solver step size of 1e-8

cset = getActiveConfigSet(hsys);
set_param(cset,'StopTime','0.015');
set_param(cset,'Solver','FixedStepAuto');
set_param(cset,'FixedStep','1e-8');
out = sim('lqitest',cset);

figure
plot(out.y),grid

Get the open loop system at the input to the plant. linearize automagically deals with the hybrid model

sysolsimulink = -1*linearize('lqitest',linio('lqitest/control',1,'looptransfer'));

Build the system using connect after explicitly discretizing the plant using the ZOH assumption.

sysp = ss(plant.A,plant.B,plant.C,plant.D,'InputName','u','OutputName',{'x1' 'x2' 'y'});
%int = tf(K_lqi(3),[1 0],'InputName','e','OutputName','ui');
Kfb = ss(K_lqi(1:2),'InputName',{'x1' 'x2'},'OutputName','ufb');
esum = sumblk('e = r - y');
usum = sumblk('u = -ui - ufb');

syspdiscrete = c2d(sysp,Ts,'zoh');
intdiscrete  = tf(Ts*K_lqi(3),[1 -1],Ts,'InputName','e','OutputName','ui');
syscdiscrete = connect(syspdiscrete,intdiscrete,Kfb,esum,usum,'r','y','u');

Compare the Bode plots of the open loop system derived by linearize() and derived by connect() followed by getLoopTransfer

figure
bode(sysolsimulink,getLoopTransfer(syscdiscrete,'u',-1),{10 1e7});

The results match.