Discretize State Space feedback controller using c2d()

I'm looking to discretize a State Space feedback controller using c2d(), to find loop rate limits with which this feedback system remains stable.

I tried what I thought was correct (below), but the discrete sys step response to approach 0 as the Ts decreases (rate increase). Also, is there a way in Matlab to create a discretized feedback controller around a continuous plant?

Background:

If I just do c2d() on the continuous CL (closed-loop) sys, as expected the discrete sys's step()/etc plots are always stable: just a sampled version of the continous CL response output. But I'm looking for controller stability, and just c2d(closed-loop sys) seems like it wouldn't show stability in presence of disturbances, say.

Is the discretization process below correct, as far as Matlab command usage?

a) Find continuous-domain plant A, B, C, D, and feedback gains K, eg:

sys = ss(A, B, C, D)
K = lqr(sys, Q, R)
sysCL = ss(sys.A - sys.B * K, sys.B * K * [1;0], sys.C, sys.D)

And to check response to disturbance inputs:

F = [0; 1/m] % since this is a 2nd-O spring/mass/damper system
sysCL = ss(sys.A - sys.B * K, [ sys.B * K * [1;0], F], sys.C, sys.D)

Use the cts version as inputs to discretization process, and as comparison

b) To find the equivalent discrete version, find discrete gains K_disc, eg:

sys_d = c2d(sys, Ts)
K_d = lqrd(sys_d.A, sys_d.B, Q, R, Ts)
sysCL_d = ss(sys_d.A - sys_d.B * K_d, sys_d.B * K_d * [1;0], sys_d.C, sys_d.D, Ts)
% Ref gain is small when loop rate is low...why do I need to scale sysCL_d by 1/dcgain(sysCL_d)?
dcgain(sysCL_d) 

And to check response to disturbance inputs:

% Is F the same for a discretized system? sysCL_d.B has elements for both states, unlike sysCL.B
F = [0; 1/m] %since this is a 2nd-O spring/mass/damper system 
sysCL_d = ss(sys_d.A - sys_d.B * K_d, [ sys_d.B * K_d * [1;0], F ], sys_d.C, sys_d.D, Ts)
% Ref gain is small when loop rate is low...why do I need to scale sysCL_d by 1/dcgain(sysCL_d)?
dcgain(sysCL_d) 

Questions:

Then I check step(sysCL, sysCL_d)

1) The above causes the discrete sys step response to approach 0 as the Ts decreases (rate increase). The cts version steps to 1 as expected.

What am I missing?

2) Is there a way in Matlab to create a discretized feedback controller around a continuous plant?

The disturbance outputs for the discrete version approach 0 as the rate goes down (Ts increase) -- but I would have expected a lower rate to cause the full dist to be at the output y until the next controller tick. Is this because only a sampled y is being displayed. I'm looking for a way to display the cts system states.

Answer 1

[Paul]

Hi John,

lqrd returns the gains for a discrete time controller. The Control System Toolbox cannot model hybrid systems, as far as I know (you can do that in Simulink if you want). With a continuous plant and discrete time controller the typical approach is to discretize the plant and then apply the feedback to approximate the hybrid, closed loop system. In other words, we can't close the loop using the continuous plant and the dicrete controller as done in the code in the Question

Here's an example.

Continuous plant

zeta = 0.3;
wn_rPs = 250;
ampl = 1;
m = 1e-5;
b = zeta * 2 * wn_rPs * m;
k = wn_rPs^2 * m;
A = [0 1; -k/m -b/m];
B = [0; 1/m];
C = [1 0];
D = 0;
plant = ss(A,B,C,D);

Desing LQR feedback gains

Q = [1000 0 ; 0 0.001];
R = 0.01;
K = lqr(plant, Q, R);

Continuous, closed loop system

syscl = ss(plant.A - plant.B*K,plant.B*K*[1;0],plant.C,plant.D);

Find the eigenvalues of the closed loop system to define an appropriate sampling period

format short e
e = eig(syscl)

e = 2×1

1.0e+00 * -1.0006e+03 -3.1605e+04

Ts = -1/e(2)/10;

Design discrete time gains

Kd = lqrd(plant.A,plant.B,Q,R,Ts);

Discretize the plant

plantd = c2d(plant,Ts);

Discrete time, closed loop system

syscld = ss(plantd.A - plantd.B*Kd,plantd.B*Kd*[1;0],plantd.C,plantd.D,Ts);

Compare step responses, they are nearly identical.

step(syscl,syscld)

[John]

Thanks @Paul. Hmmm, that's the same process, then, that I had.

If the rates are lowered (eg going all the way down to an unreasonable 100 Hz), in your code you'd also see the step response steady-state --> 0, no? This doesn't make sense to me.

Similarly, disturbance response acts differently, even with a high update rate system -- even when the two ref responses approach the same value, as you'd shown.

Perhaps there disturbance vector needs to be changed for the discrete system somehow? Per this comment:

% Is F the same for a discretized system? sysCL_d.B has elements for both states, unlike sysCL.B
F = [0; 1/m] %since this is a 2nd-O spring/mass/damper system 

Thanks -- good to know about the hybrid system needing SL.

[Paul]

Yes, that is the same process you showed (for the control, but not for the disturbance). Not sure why I thought otherwise.

However, the plant has to include the disturbance input when it's discretized, which is not what you had.

zeta = 0.3;
wn_rPs = 250;
ampl = 1;
m = 1e-5;
b = zeta * 2 * wn_rPs * m;
k = wn_rPs^2 * m;
A = [0 1; -k/m -b/m];
B = [0; 1/m];
C = [1 0];
F = [0; 1/m];
D = 0;

Include the disturbance input here

plant = ss(A,[B F],C,D);

Q = [1000 0 ; 0 0.001];
R = 0.01;
K = lqr(plant(1,1), Q, R);
syscl = ss(plant.A - plant.B(:,1)*K,[plant.B(:,1)*K*[1;0] plant.B(:,2)],plant.C,plant.D);
e = eig(syscl)

e = 2×1

1.0e+04 * -0.1001 -3.1605

Ts = -1/e(2)/10;

Kd = lqrd(plant.A,plant.B(:,1),Q,R,Ts);

Discretize plant including disturbance and then close the loop

plantd = c2d(plant,Ts);
syscld = ss(plantd.A - plantd.B(:,1)*Kd,[plantd.B(:,1)*Kd*[1;0] plantd.B(:,2)],plantd.C,plantd.D,Ts);
step(syscl(1,2),syscld(1,2))

You'd see the responses to the disturbance would match better with Ts smaller.

I'm not sure what the concern is with making Ts larger. As Ts gets larger, the discretized plant model becomes a worse approximation to the continuous plant and the Kd gains are designed for a discrete time plant model that doesn't really approximate the continuous plant, so there's no reason to expect the closed loop response of the discretized system to approximate the closed loop response of the continuous system if Ts is too large.

[John]

@Paul

plant = ss(A,[B F],C,D);
...
Kd = lqrd(plant.A,plant.B(:,1),Q,R,Ts);

Ah! Thanks. That makes sense; the disturbance needs to be variable size with variable step, and that can only be accomplished as part of a full system needed as the input for discretization.

By the way, I think the two inputs of B can be referenced in this way as well instead of plant.B(: , 1), perhaps a little more clearly:

plant(1).B for u

plant(2).B for F

"I'm not sure what the concern is with making Ts larger. As Ts gets larger, the discretized plant model becomes a worse approximation to the continuous plant and the Kd gains are designed for a discrete time plant model that doesn't really approximate the continuous plant, so there's no reason to expect the closed loop response of the discretized system to approximate the closed loop response of the continuous system if Ts is too large."

I see. I assumed that as the approximation is worse, that doesn't matter that it doesn't approximate the real plant. I would have thought what matters is that the found gains bring that plant to 1 with a ref step. So wouldn't the gain computation find whatever gains are needed to bring that plant (even if it's not an accurate representation) to the correct steady-state value when feedback is closed? The gains are found with the discrete system via lqrd(), so it already knows how that plant behaves...

[Paul]

Did you try

plant(1).B for u

plant(2).B for F

I don't think that will work.

Even in the LQR continuous case, the found gains don't bring the ouput to 1 with a unit step input because there is no free integrator in the open loop transfer function at the error signal.

zeta = 0.3;
wn_rPs = 250;
ampl = 1;
m = 1e-5;
b = zeta * 2 * wn_rPs * m;
k = wn_rPs^2 * m;
A = [0 1; -k/m -b/m];
B = [0; 1/m];
C = [1 0];
D = 0;
plant = ss(A,B,C,D);
Q = [1000 0 ; 0 0.001];
R = 0.01;
K = lqr(plant, Q, R);
syscl = ss(plant.A - plant.B*K,plant.B*K*[1;0],plant.C,plant.D);
dcgain(syscl)

ans = 0.9980

So not exactly unity. Analytically, the dc gain is

K(1)/(K(1) + k)

ans = 0.9980

So it would be unity if not for the spring constant.

If the spring constant is 0, then the plant has a pole at s = 0 and the discretized plant model will have a pole at z = 1, and we'd get the unity steady state tracking in both continous and discrete time.

plant.A(2,1) = 0;
syscl = ss(plant.A - plant.B*K,plant.B*K*[1;0],plant.C,plant.D);
Ts = 0.01;
Kd = lqrd(plant.A,plant.B,Q,R,Ts);
plantd = c2d(plant,Ts);
syscld = ss(plantd.A - plantd.B*Kd,plantd.B*Kd*[1;0],plantd.C,plantd.D,Ts);

step(syscl,syscld)

[John]

@Paul

1) "Did you try plant(1).B for u, plant(2).B for F

I don't think that will work."

Hmmm, I had thought I worked but could be wrong... In Matlab, isn't a two-input SS system a vector of two SS systems?

Here's my output. Note the updated F to be 2*B, just for debugging clarity (else columns of u would be the same ).

A = [0 1; -k/m -b/m];
B = [0; 1/m];
C = [1 0];
F = [0; 2*1/m];
D = 0;
plant = ss(A,[B F],C,D)
plant(1).B
plant(2).B

The output is

plant =
 
  A = 
              x1         x2
   x1          0          1
   x2  -6.25e+04       -150
 
  B = 
          u1     u2
   x1      0      0
   x2  1e+05  2e+05
 
  C = 
       x1  x2
   y1   1   0
 
  D = 
       u1  u2
   y1   0   0
 
Continuous-time state-space model.

ans =

0
1.0000e+05


ans =

0
2.0000e+05

which seems like it's working.

2) "So it would be unity if not for the spring constant.

If the spring constant is 0, then the plant has a pole at s = 0 and the discretized plant model will have a pole at z = 1, and we'd get the unity steady state tracking in both continous and discrete time."

I see -- thanks. That's starting to make sense, but I'll have to think about it a bit more. If I can apply a manual gain to lqrd()'s results, and that gain can bring the response to step to 1 (same for lqr()), I wonder why the A.R.E can't output gains to step nearer 1 with lower sampling rates (I know it's a specific computation with certain assumptions about the system, so perhaps those assumptions break down as Ts approaches the time constant of the system)

3) Related to the original question, when I try your approach: I do get the same result as you -- thanks for helping with that.

If I were to go through the same process with a more complex continuous model, for example LQRI (which is working in cts-domain), the process would be the same as above for all diff equations, right?

That is, say I find discrete gains Kd = [Kx_lqgi_d, Ki_lqgi_d], and discrete observer gains L_lqgi_d. As before, can I construct A_closedLoop_discrete like the following below, with the equations being the direct analogue to the cts version?

Ignore whether the equations themselves are right; I'm asking just about the principal that plant --> plantd, and gains --> gainsd, etc, regardless of how complex the A_closed_loop construction is.

I'm asking this since when c2d() maps of A --> Ad, it includes things like identity matrices and Ts applied to A, so I'm unsure if I can construct a larger A per below -- but i think this is correct because each differential equation (state level) is already made of discretized components.

A11 = plantd.A % cts is: plant.A
A12 = zeros(2, 1) % cts is: zeros(2, 1)
A13 = -plantd(1).B * Kx_lqgi_d % cts is: -plant(1).B * Kx_lqgi
A14 = -plantd(1).B * Ki_lqgi_d % cts is: plant(1).B * Ki_lqgi
A31 = L_lqgi_d * plantd.C % cts is: L_lqgi * plant.C 
...
A33 = plantd.A - plantd(1).B * Kx_lqgi_d - L_lqgi_d * plantd.C % etc...
...
A43 = -plantd.C;
A44 = 0;
A_closedLoop_d = [    A11, A12, A13, A14;
                      A21, A22, A23, A24;
                      A31, A32, A33, A34;
                      A41, A42, A43, A44   ];
...
sysCl_d = ss(A_closedLoop_d, ...
    [plantd(1).B * Kx_lqgi_d* [1;0] plantd(2).B], ...
    plantd.C, plantd.D);
% cts is:
% sysCl = ss(A_closedLoop, ...
%    [plant(1).B * Kx_lqgi* [1;0] plant(2).B], ...
%    plant.C, plant.D);

[Paul]

1). Looks like it does work. Learn something new every day.

2). I'm not exactly sure what you're asking here. lqr or lqrd return the optimal gains based on the input parameters Q, R, and Ts if applicable. You can always to the dcgain adjustment outside the loop after the fact if desired.

3). I'm not following those equations. If you design the control gains Kd and the observer gains Ld both based on the same, discretized plant model, plant_d, then the closed loop system should be constructed in the usual way based on plant_d, the discrete observer equations, and the control.

[John]

Thanks @Paul. Makes sense for (2).

For (3), let me try asking from a different angle since i wasn't very clear...

In the LQR example, the augmented ([ref F]) plant system was discretized first. This made sense because this is a physical system representation.

Let's say I'm next trying to do LQI discretization. In this case the reference comes through the integrater, which is not a physical plant state. (state vector x: [x1, x2, xi]. Then u = [0; 0; 1] for a step input).

Normally, for the continuous case, it would be:

% x = [x1 (pos), x2 (vel), xi]
% Reference input at integrator (x3), and unity gain (unlike ref inputs to
% vel level that are not unity gain)
B_lqi = [0; 0; 1]
C_lqi = [physSys.C, 0];
uLqi = [ref distVel];
% Dist at vel state level, ie acceleration
Tau_lqi = [0; 1 / J_kgm2_phy; 0];
sysCl_Lqi = ss(Acl_lqi, [B_lqi, Tau_lqi], C_lqi, 0);
sysCl_Lqi.StateName = {'theta', 'angVel', 'xi'};
sysCl_Lqi.InputName = {'r(ref)', 'd(F)'};
sysCl_Lqi.OutputName = {'theta'};

But it seems like in the discretized workflow above in your post, additional states must be included in an augmented plant before the closed loop system is built witrh ss(). In other words, a discretized B in ss() needs to already be expressable as [sysd(1).B, sysd(2).B], specifically with cts sys.B coming from an augmented plant. Eg for LQI, LQGI, etc, it seems like a plant needs be built so that B includes all controller states, which seems counterintuitive:

How would this generalize to something with arbitrary states whose structure is normally added after the plant is built, since the designer has specific knowledge of the control structure, eg xi or observer states? The plant is just the physical system itself without knowledge of the controller...

How would a discretized LQI be created in that c2d() workflow, eg using the code framework above?

[Paul]

If doing direct, discrete time integral control design, IMO it doesn't make sense to augment the plant with a continuous integrator and then do c2d on the augmented plant. Instead, I would do c2d on the physical plant to get a discrete time physical plant model, and then augment the discrete time plant model with a discrete integrator and go from there to do the design. For implementation, the integrator in the compensation should be of the same form as that used in the design (e.g., Euler, Tustin, whatever).

I think that's why the doc page for lqi specifies the form of the discrete-time integrator, so that the user can implement in the form that it was assumed for design of the gains.

Having said that, it wouldn't suprise me if using c2d on the augmented plant and then designing the gains can be made to work in some sense, but you'd have to figure out how to relate the ZOH approximation used in c2d (by default) to how the integral control would actually be implemented.

[John]

@Paul Thanks. "Instead, I would do c2d on the physical plant to get a discrete time physical plant model, and then augment the discrete time plant model with a discrete integrator and go from there to do the design. For implementation, the integrator in the compensation should be of the same form as that used in the design (e.g., Euler, Tustin, whatever)."

I see, thanks; good point on integrator consistency. But, when I formulate the system in SS, it seems there isn't an integrator -- it's a derivative with a straight subtraction: xi_dot = ref - y. Where do you mean to add that, in this context?

Augmented plant per your method:

% LQI states: x = [x1, x2, xi]
F = [0; 1 / J];
% Plants include dist: B --> [B, F]
plant = ss(A, [B, F], C, D)
plantd = c2d(plant, Ts)

Discrete LQI:

% Same Q, R as cts version
[Kd, S1, P1] = lqi(plantd(1), Q, R)
% [A-BKx, -BKi;
%   -C,   0]
Ad_cl = [plantd.A - plantd(1).B * Kd(1:2), -plantd(1).B * Kd(3);
                    -plantd.C,                        0]

% Input to integrator level
Bd = [0; 0; 1];
Cd = [plantd.C 0];
Dd = [0];
% Augment F to add 0 at xi level. Rest is just plantd's dist part of B
F_aug_d = [plantd(2).B; 0]
sysCl_d = ss(Ad_cl, [Bd F_aug_d], Cd, Dd, Ts);

Cts LQI version:

% Same Q, R as disc version
[K, S1, P1] = lqi(plant(1), Q, R)
% [A-BKx, -BKi;
%   -C,   0]
Acl = [plant.A - plant(1).B * K(1:2), -plant(1).B * K(3);
    -plant.C,                         0]

% Input to integrator level
B = [0; 0; 1];
C = [plant.C 0];
D = [0];
% Augment F to add 0 at xi level. Rest is just plantd's dist part of B
F_aug = [plant(2).B; 0]
sysCl = ss(Acl, [B F_aug], C, D, Ts);

This doesn't give the same Ref step() response as the discrete version:

Zoomed in time only:

Settling at 1 seems correct, but x100 faster response. I also tried scaling the integrator (xi) level by Ts, in case that's what you meant, but couldn't find a way to get them to match.

It's unclear to me: how would I use Matlab to discretize basic state-space systems like LQI, for systems where the continuous version is working well?

[Paul]

Here's an example. I'm using connect here, because I still think it's easier than doing all matrix manipulations.

Plant model

overShootPec = 5;
zeta = abs(log(overShootPec/100)) / sqrt(pi^2 * log(overShootPec/100)^2);
wn_rPs = 2 * pi * 40;
%% OL system
ampl = 1;
J_kgm2 = 1e-5;
b = zeta * 2 * wn_rPs * J_kgm2;
k = wn_rPs^2 * J_kgm2;
A = [0 1; -k/J_kgm2 -b/J_kgm2];
B = [0; 1/J_kgm2];
C = [1 0];
D = [0];
massSys = ss(A, B, C, D, ...
    'StateName', {'theta', 'theta_dot'}, ...
    'InputName', {'u'}, ...
    'OutputName', {'theta'});
%Plant with disturbance
% define F so that disturbance is a torque
F = [0; 1/J_kgm2];
sysOlWithDist = ss(massSys.A, [massSys.B  F] , eye(2),0);
sysOlWithDist.StateName = {'theta', 'theta_dot'};
sysOlWithDist.InputName = {'u', 'd'};
sysOlWithDist.OutputName = {'theta','theta_dot'};

Q and R for LQR design

H = [2*0.7*1000 1 -1000^2];
Q = H.'*H;
R = 100;

Design with integral control using lqi.

Kc = lqi(sysOlWithDist(1,1), Q, R);

Continuous, closed-loop system

sysintcont = tf(1,[1 0],'OutputName','xi','InputName','e');
K = ss(-Kc,'InputName',{'theta','theta_dot','xi'},'OutputName','u');
S = sumblk('e = r - theta');
sysc = connect(sysOlWithDist,K,S,sysintcont,{'r','d'},'theta');

Discrete time integrator, using Tustin's rule

Ts = 1e-5;
sysintdisc = c2d(sysintcont,Ts,'Tustin');

Discrete time, closed loop system by discretizing the plant, and connecting with the discrete integrator, still using the Kc gains.

sysd1 = connect(c2d(sysOlWithDist,Ts,'zoh'),K,S,sysintdisc,{'r','d'},'theta');

Now, do direct discrete design with integral control.

First, define the integrator. Here, I'm using the same integrator implementation as in lqi. Rethinking my comment from above, it will be more difficult, if possible, to use a discrete time integrator that has direct feedthrough, Also, I'm adding a negative sign because I want to use the same Q as used above in the call to lqi, and as we've previously discussed lqi puts adds a negative sign because of the negative feedback it assumes for the error signal

sysintdisc = ss(-tf(Ts,[1 -1],Ts));

Ensure that state variable is the same as the ouptut, which is the (negative of the) integral of the input.

sysintdisc.B = sysintdisc.B*sysintdisc.C;
sysintdisc.C = 1;
set(sysintdisc,'InputName','theta','OutputName','xi','StateName','xi');

Define the discrete time, augmented plant

augplant = connect(c2d(sysOlWithDist(1,1),Ts,'zoh'),sysintdisc,'u','xi');

Compute the gains.

Kd = lqr(augplant,Q,R);

Redefine the intgegrator for implementation; it has the same transfer function as for the design, except for the negative sign that is now captured in the formation of the error.

sysintdisc = ss(tf(Ts,[1 -1],Ts));
set(sysintdisc,'InputName','e','OutputName','xi','StateName','xi');
K = ss(-Kd,'InputName',{'theta','theta_dot','xi'},'OutputName','u');
sysd2 = connect(c2d(sysOlWithDist,Ts,'zoh'),K,S,sysintdisc,{'r','d'},'theta');

Compare all three responses

step(sysc(1,1),sysd1(1,1),sysd2(1,1));

I think this line in your code might be problematic

% Ad_cl = [plantd.A - plantd(1).B * Kd(1:2), -plantd(1).B * Kd(3);
%                     -plantd.C,                        0]

In discrete time we have: y(k) = C*x(k) + D*u(k)

According to the lqi doc, the integrator state equation is

xi(k+1) = xi(k) + T*(r(k) - y(k))

or (assuming D = 0 as is the case in this model)

xi(k+1) = xi(k) -T*Cx(k) + T*r(k)

so the second row should be

% Ad_cl = [plantd.A - plantd(1).B * Kd(1:2), -plantd(1).B * Kd(3);
%                     -plantd.C*Ts,                        1]

And you'll need to multiply Ts wherever r(k) comes in on the inputs.

[John]

Thanks @Paul. I'll need to digest most of what you wrote, but will reply to this part in the meantime:

"According to the lqi doc, the integrator state equation is

xi(k+1) = xi(k) + T*(r(k) - y(k))

so the second row should be

% Ad_cl = [plantd.A - plantd(1).B * Kd(1:2), -plantd(1).B * Kd(3);
%                     -plantd.C*Ts,                        1]

"

Ah, thanks for catching that. That does fix it.

Good point: since all xn(k+1) = xn(k) + ... , the entire A matrix needs +eye(3), which was already added to the first two levels via plant.A --> plantd.A, but wasn't added to the last row since the 3rd row is built manually. Similarly, any diff eq from cts --> disc needs a Ts multiplier when manually doing the conversion (coming from whichever discrete derivative estimation technique is used), so *Ts needed (as you said) on elements that row as well. Thanks!

I'll use connect() as verification; but since I'm building discrete models in Simulink, I need to understand how to manually build them up using Matlab.

[John]

@Paul There's something unexpected with the formulation above. The step responses (r) are as expected (as you showed), but the dist resp (F) are not.

Specifically, as discretization (loop) rates increase, the discretized version unexpectedly seems to become more effective than the cts version. The discretized Dist response doesn't converge to Cts, but passes through that response and becomes progressively smaller. Here are three snapshots.

As expected, disc response worse at a low rate ~10 kHz.

Now apparently disc is more effective than cts at 100 kHz

And disc continues to improve as rate increases:

Of course, this is some kind of simulation artifact, since continuous-time (inf loop rate) is the best-case optimal scenario that can't be exceeded.

c2d() is using tustin, and other methods don't change this result.

Is this a Matlab limitation, eg rounding error accumulation? If so, 100 kHz seems slow for that to occur, so it would be surprising.

Or how is this possible? Step() should automatically take steps according to the discrete model's Ts that's inputted.

[Paul]

Simulink can (should?) be used as graphical representation of connect. Personally, I think it would be much better to use Simulink that way, i.e., have different blocks for the plant, control, and observer. That way you model a hybrid system (continous plant, discrete control and observer), can log inputs/outputs of each element in the loop, add io analysis points if needed for linear analysis, etc. I think you lose a lot of visbility into the system by trying to manually create a single, state space model for the entire closed loop system. Of course, your mileage may vary.

[John]

100% agreed; that's my process as well.

I'll clarify: the step after doing what you described (connecting blocks of varying rates, diagnostics, etc) is to eventually implement this on a physical system. Given what that system is, this process will be manually recreating the simulink / matlab system, so understanding what ML / SL does is a needed step... Using ML/SL to create a working system at a high level, is certainly part of that process.

Hence being a bit unclear on what Matlab simulating, if increasing loop rates lead to performance better than a continuous sytems (which shouldn't be possible).

[Paul]

I'm seeing basically the same response to a step disturbance for all three models as long as Ts is small enough. Same code as above, just added a plot for the repsonse to the disturbance and cycled through three values of Ts

overShootPec = 5;
zeta = abs(log(overShootPec/100)) / sqrt(pi^2 * log(overShootPec/100)^2);
wn_rPs = 2 * pi * 40;
%% OL system
ampl = 1;
J_kgm2 = 1e-5;
b = zeta * 2 * wn_rPs * J_kgm2;
k = wn_rPs^2 * J_kgm2;
A = [0 1; -k/J_kgm2 -b/J_kgm2];
B = [0; 1/J_kgm2];
C = [1 0];
D = [0];
massSys = ss(A, B, C, D, ...
    'StateName', {'theta', 'theta_dot'}, ...
    'InputName', {'u'}, ...
    'OutputName', {'theta'});
%Plant with disturbance
% define F so that disturbance is a torque
F = [0; 1/J_kgm2];
sysOlWithDist = ss(massSys.A, [massSys.B  F] , eye(2),0);
sysOlWithDist.StateName = {'theta', 'theta_dot'};
sysOlWithDist.InputName = {'u', 'd'};
sysOlWithDist.OutputName = {'theta','theta_dot'};

H = [2*0.7*1000 1 -1000^2];
Q = H.'*H;
R = 100;
Kc = lqi(sysOlWithDist(1,1), Q, R);

sysintcont = tf(1,[1 0],'OutputName','xi','InputName','e');
Kcc = ss(-Kc,'InputName',{'theta','theta_dot','xi'},'OutputName','u');
S = sumblk('e = r - theta');
sysc = connect(sysOlWithDist,Kcc,S,sysintcont,{'r','d'},'theta');

for Ts = [1e-4 1e-5 1e-6]
    
sysintdisc = c2d(sysintcont,Ts,'Tustin');

sysd1 = connect(c2d(sysOlWithDist,Ts,'zoh'),Kcc,S,sysintdisc,{'r','d'},'theta');

sysintdisc = ss(-tf(Ts,[1 -1],Ts));
sysintdisc.B = sysintdisc.B*sysintdisc.C;
sysintdisc.C = 1;
set(sysintdisc,'InputName','theta','OutputName','xi','StateName','xi');

augplant = connect(c2d(sysOlWithDist(1,1),Ts,'zoh'),sysintdisc,'u','xi');
Kd = lqr(augplant,Q,R);

sysintdisc = ss(tf(Ts,[1 -1],Ts));
set(sysintdisc,'InputName','e','OutputName','xi','StateName','xi');
Kcd = ss(-Kd,'InputName',{'theta','theta_dot','xi'},'OutputName','u');
sysd2 = connect(c2d(sysOlWithDist,Ts,'zoh'),Kcd,S,sysintdisc,{'r','d'},'theta');

figure
step(sysc(1,1),sysd1(1,1),sysd2(1,1));
text(0,1,"Ts = " + string(Ts))

figure
step(sysc(1,2),sysd1(1,2),sysd2(1,2));
text(0,0,"Ts = " + string(Ts))

end