lsim() vs step() : are different responses expected?

Are different responses expected if running lsim() and step(), both with the same unit step inputs?

I have a state-space system representing a closed-loop system, that gives the expected output with a unit step input to lsim() -- but not with a unit step input to step().

Even if the system itself is incorrectly built, i'd have expected equivalent responses from step() and lsim().

Background:

The system is attached in sys.mat, and is shown at the end of the post as well. The .m file that plots these is also attached.

The states are: x = [x1, x2, xi, x1^, x2^, xi^], where ^ are estimated states. x1,x2 are pos/vel, and xi is the integral action state.

x1, x1^ , and x2, x2^ , and xi^ are all as expected.

But xi is not as expected. For some reason, with the same ref step input of 1,

step() gives an unexpected xi=0, and

lsim() gives the expected xi .

Questions:

What's the difference between step() and lsim() -- if the same system and same ref inputs are used -- that would cause outputs to differ?

Are step() and lsim() both inputting the ref to u? If they are both setting u = 1, then shouldn't sys.B provide inputs to both xi and xi^ level even when using step(), since B is [0;0;1, 0;0;1]?

For step(): see 2nd subplot, blue line = 0.

For lsim(): both plots seem as expected.

The details are perhaps not important: regardless of whether the system is correctly built or not, i'd have expected equivalent responses from step() and lsim()

But in case it matters, in the state space sys, the two relevant rows are xi, and xi^:

xi_dot = -C*x + 1*u

xi^_dot = -C*x^ + 1*u

with u = reference step input, and -C*x = y (theta) and -C*x^ = y^ (theta^)

Here's the system. It's also in the attached .mat file.

sysCl_lqgi =

A = 
              theta      angVel          xi      theta^     angVel^         xi^
theta             0           1           0           0           0           0
angVel   -6.317e+04        -160           0  -4.468e+07  -6.376e+04   1.596e+10
xi               -1           0           0           0           0           0
theta^    4.212e+04           0           0  -4.212e+04           1           0
angVel^   4.398e+08           0           0  -4.845e+08  -6.392e+04   1.596e+10
xi^               0           0           0          -1           0           0

B = 
             u      F
theta        0      0
angVel       0  1e+05
xi           1      0
theta^       0      0
angVel^      0      0
xi^          1      0

C = 
         theta   angVel       xi   theta^  angVel^      xi^
theta        1        0        0        0        0        0

D = 
       u  F
theta  0  0

Continuous-time state-space model.

Answer 1

[Paul]

@John

I was able to recreate your strange result on R2022a.

After some investigation, I think I know what's happening.

load sys.mat
sys = sysCl_lqgi;

The six states of this system are:

sys.StateName.'

ans = 1×6 cell array
    {'theta'}    {'angVel'}    {'xi'}    {'theta^'}    {'angVel^'}    {'xi^'}

After examining the state space matrices and based on the state names, it looks like you started with a 2-state model, augmented it with an integrator, then designed an observer for the augmented plant (including the integrator) and then closed the loop through the control applied to the estimated states, where the external reference signal forms the error that feeds xidot and xihatdot. However, the third column of sys.A is all zeros, so even if xi is excited by the reference input, xi will have no effect on the output, which is theta. Same thing with xi being excited by the disturbance input. Consequently, xi can be eliminated from the system with no impact on the input/output response. This elimination is called structural minimal realization sminreal.

msys = sminreal(sys);
msys.StateName.'

ans = 1×5 cell array
    {'theta'}    {'angVel'}    {'theta^'}    {'angVel^'}    {'xi^'}

Note that sminreal eliminates the third state.

I'm 99% certain (can't be 100% because I ran into a function that was implemented in p-code and so I couldn't examine it) that deep in the bowels of the code that computes the step response, it does a structural minimal realization (after converting the input sys to its discrete time approximation) of the system. But, in the outputs returned to the user, the eliminated state(s) is filled with zeros, so as to keep the dimensions of the state vector output of step() aligned with the definitions of sys. I don't have 2023a installed locally and so can't see what changed in 2023a.

You can test this by changing the C matrix to

sys.c = [0 0 1 0 0 0];

In so doing, that third state can't be eliminated because it's needed for the output

msys = sminreal(sys);
msys.StateName.'

ans = 1×6 cell array
    {'theta'}    {'angVel'}    {'xi'}    {'theta^'}    {'angVel^'}    {'xi^'}

You'll see that if you run your code with this sys that you get the step response plots you expect for the state vector.

[John]

Interesting, @Paul. Thanks for looking more into it.

"However, the third column of sys.A is all zeros, so even if xi is excited by the reference input, xi will have no effect on the output, which is theta. Same thing with xi being excited by the disturbance input. Consequently, xi can be eliminated from the system with no impact on the input/output response."

Yes, exactly. This is by design; that state is unused, but I have a few formulation variations that switch back and forth between using xi and xi^ -- so that specific implementation doesn't use it.

"But, in the outputs returned to the user, the eliminated state(s) is filled with zeros, so as to keep the dimensions of the state vector output of step() aligned with the definitions of sys. I don't have 2023a installed locally and so can't see what changed in 2023a."

Huh, that seems like a good thought. Agreed, then, that something seems like something changed from 2022b and 2023a. But, in that case, niavely it still seems like the output should not be 0s, since the SS equations still exist at that level, even though it doesn't impact the output. Overall, that seems plausabile (if confusing why that choice would have been made)

I wonder what Mathworks' investigation will show...

[Paul]

My initial reaction was the same as yours, i.e., the full state vector should be returned to the user even if one or more states can be structurally eliminated. Curious to know why that design decision was originally made. Maybe there was some corner case where keeping all of the states had some impact on the numerical accuracy of the solution. Though, in that regard, I found it curious that the state was eliminated only after the model was converted to discrete time. Anyway, assuming my assertion is correct, my guess is we'll never find out why, only that it's been addressed in 2023a and the workaround for prior versions is to ensure that sys is structurally minimal.

[John]

@Paul @Sam Chak @Walter Roberson Mathworks' ticket reply is the same as what Paul found. In short, 2023a doesn't have this issue (as suspected), and the accidental removal is as Paul described above. Giving the answer here for anyone else who runs across this:

"Thank you for your patience in awaiting my reply. The developers would like to sincerely thank you for bringing this long-standing bug to our attention.

After loading "sys.mat" into the workspace, please enter "sysCl_lqgi" into the command line. You will observe that the third column of A contains only 0's, and the third column of C is 0. This means that "xi" affects no other states and does not affect the output. Thus, it can be eliminated when calculating the output. This also manifests itself in a pole-zero pair in the s-domain. Afte running "zpk(sysCl_lqgi)", observe the pole/zero pair in the transfer function from "u" to "theta".

Prior to R2023a, the "step(sys)" algorithm calculated the output based on the system returned from "sminreal(sys)". Without going into details, I can say that the new algorithm does not perform this simplification.

As a workaround in R2022b, you can insert the commands:

>>> [a,b,c,d] = ssdata(sys);

>>> sys = ss(a,b,eye(6),0);

before the command:

>>> [y,tout,x] = step(sys);

This will add extra outputs, each dependent on a different state. Thus, "hard 0" based minimization of state-space models will not affect the output.

Please let me know if this satisfies your workflow, and whether you have any further questions."

[Paul]

I think a better approach would be

sys = ss(a,b,[c;eye(size(a))],[d;zeros(size(a,1),size(b,2))]);
[y,tout,x] = step(sys);

which preserves the original outputs specified by the user as the first elements in y.

[Paul]

I wonder why the update in 2023a is not listed as bug fix or as a change in behavior in the Release Notes

Answer 2

[Paul]

The code and the system in the .mat file does not recreate the step plot. There is also an undefined variable.

load sys.mat
% step()
sys = sysCl_lqgi;
[y,tout,x] = step(sys);

figure
subplot(2, 1, 1);
plot(tout, x(:, 1, 1), 'b', 'LineWidth', 1.5);
hold on
plot(tout, x(:, 4, 1), 'k:', 'LineWidth', 2.5)
grid on;
title('x1', 'FontSize', 14)
legend({'$x_1$', '$\hat{x_1}$'},'Interpreter','latex', 'FontSize', 16);
xlim([0, 0.02]);

subplot(2, 1, 2);
plot(tout, x(:, 3, 1), 'b', 'LineWidth', 1.5);
hold on
plot(tout, x(:, 6, 1), 'k:', 'LineWidth', 2.5);
grid on;
title('xi', 'FontSize', 14)
legend({'$x_i$', '$\hat{x}_i$'},'Interpreter','latex', 'FontSize', 16);
xlim([0, 0.02]);
xlabel('Time [s]', 'FontSize', 14)
sgtitle('STEP(): Ref response');

% lsim()

step inputs a step simultaneously to both inputs. So the second column of u should also be ones, not zeros.

u = [ones(size(t)) zeros(size(t))];

Unrecognized function or variable 't'.

[y, tout, x] = lsim(sys, u, t, zeros(6,1));

figure;
subplot(2, 1, 1);
plot(t, x(:, 1, 1), 'b', 'LineWidth', 1.5);
hold on;
plot(t, x(:, 4, 1), 'k:', 'LineWidth', 2.5);
grid on;
title('x1')
legend({'$x_1$', '$\hat{x_1}$'},'Interpreter','latex', 'FontSize', 16);
xlim([0, 0.02]);

subplot(2, 1, 2);
plot(t, x(:, 3, 1), 'b', 'LineWidth', 1.5);
hold on;
plot(t, x(:, 6, 1), 'k:', 'LineWidth', 2.5);
grid on;
title('xi', 'FontSize', 14);
legend({'$x_i$', '$\hat{x}_i$'},'Interpreter','latex', 'FontSize', 16);
xlim([0, 0.02]);
xlabel('Time [s]', 'FontSize', 14)
sgtitle('LSIM(): Ref response', 'FontSize', 14)

[John]

Hi @Paul. That's odd. I run the same code, and do not get the plot you show where xi tracks xi^ (ie nonzero xi). Mine still seems to show xi = 0.

I cleared all workspace and closed all figures.

I made a new file and copied in your code. I used a downloaded .mat file that was attached to your post (to make sure nothing changed from what you used), and I still get the first step() plot I showed where xi = 0.

I'm attaching the .mat and .m file I used -- should be same as yours.

What happens when you run what's attached to this post?

[Paul]

Here you go. I just copied the code out of testPaulStepPost.m and ran it here with no changes.

clear all
close all

load sys.mat
% step()
sys = sysCl_lqgi;
[y,tout,x] = step(sys);

figure
subplot(2, 1, 1);
plot(tout, x(:, 1, 1), 'b', 'LineWidth', 1.5);
hold on
plot(tout, x(:, 4, 1), 'k:', 'LineWidth', 2.5)
grid on;
title('x1', 'FontSize', 14)
legend({'$x_1$', '$\hat{x_1}$'},'Interpreter','latex', 'FontSize', 16);
xlim([0, 0.02]);

subplot(2, 1, 2);
plot(tout, x(:, 3, 1), 'b', 'LineWidth', 1.5);
hold on
plot(tout, x(:, 6, 1), 'k:', 'LineWidth', 2.5);
grid on;
title('xi', 'FontSize', 14)
legend({'$x_i$', '$\hat{x}_i$'},'Interpreter','latex', 'FontSize', 16);
xlim([0, 0.02]);
xlabel('Time [s]', 'FontSize', 14)
sgtitle('STEP(): Ref response');

[John]

Thanks @Paul. Okay, it seems something is perhaps not right with Matlab: I don't get the same output with the same code and variable.

I tried Impulse() as well, with the same result of xi = 0.

So, maybe there's something in how Matlab interprets the variable, vs the step() function having a bug.

I posted due to spending hours trying to track down where i was going wrong with using step()...thankfully it seems (perhaps) nothing on my end...

What's the next step to debug this? Should I file a service request?

If yes, as info that I can put into the service request, what platform are you running on?

(I'm OSX, running 2022b)

[Paul]

Before doing anything, what output do you get from these commands?

load sys.mat
% step()
sys = sysCl_lqgi;
which step(sys);

/MATLAB/toolbox/shared/controllib/engine/@DynamicSystem/step.m  % ss method

[John]

I run (with the addition of clear all)

clear all
load sys.mat
% step()
sys = sysCl_lqgi;
which step(sys);

and get:

/Applications/MATLAB_R2022b.app/toolbox/shared/controllib/engine/@DynamicSystem/step.m % ss method

[Paul]

I have no idea why you'd be getting such a different result. To answer your question from above, I'm running everything right here, on Answers. If you do file a service request, you can inlude a link to this thread.

[John]

Thanks @Paul. Oh I see; I just ran it here on Answers as well just now, and got the rame result as you (ie plots as expected). I'll file a service request and include a link to this thread, and update when i know more.

[Sam Chak]

Hi @John

It's highly probable that you have graphed one of the state responses from Step Input #2, , but you may not have realized it.

[John]

@Sam Chak, wouldn't @Paul have seen the same output, then, with the same code? I think the confusing thing here is that I'm getting two different results: if I Run here in Answers like Paul did then I see non-0's, or if i run on my local machine I see all 0's. And none of the other states are 0 -- only xi is.

If I'm plotting the wrong output, wouldn't that occur here in Answers as well, and for Paul as well?

Asking since I might be misunderstanding what you mean :)

"It's highly probable that you have graphed one of the state responses from Step Input #2, , but you may not have realized it."

Do you mean the F input, vs control u input? (If so, that was intentional: disturbance input)

If not, could you elaborate on what you mean by "u2"? Which variable is that in the state space formulation i posted?

[Sam Chak]

Hi @John,

Can you check the transfer function of on your machine using my code below? should the 3rd state from input 1 to output.

load('sys.mat')
A   = sysCl_lqgi.A;
B   = sysCl_lqgi.B;
C   = eye(6);
D   = 0;
sys = ss(A, B, C, D);
G   = tf(sys)

G =
 
  From input 1 to output...
                         1.596e10 s^2 + 6.747e14 s + 7.127e18
   1:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                 1.596e10 s^3 + 6.747e14 s^2 + 7.127e18 s - 5.948e-11
   2:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
          s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.044e13 s^2 + 1.998e16 s + 537.4
   3:  ----------------------------------------------------------------------------
       s^6 + 1.062e05 s^5 + 3.194e09 s^4 + 3.045e13 s^3 + 2.065e16 s^2 + 7.127e18 s
 
                         1.596e10 s^2 + 6.747e14 s + 7.127e18
   4:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                 1.596e10 s^3 + 6.747e14 s^2 + 7.127e18 s - 5.564e06
   5:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
              s^4 + 1.062e05 s^3 + 3.194e09 s^2 + 3.044e13 s + 1.998e16
   6:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
  From input 2 to output...
                    1e05 s^3 + 1.06e10 s^2 + 3.177e14 s + 1.596e15
   1:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                  1e05 s^4 + 1.06e10 s^3 + 3.177e14 s^2 + 1.596e15 s
   2:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                     -1e05 s^3 - 1.06e10 s^2 - 3.177e14 s - 1.596e15
   3:  ----------------------------------------------------------------------------
       s^6 + 1.062e05 s^5 + 3.194e09 s^4 + 3.045e13 s^3 + 2.065e16 s^2 + 7.127e18 s
 
                          4.212e09 s^2 + 3.132e14 s - 0.2129
   4:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                         4.398e13 s^2 - 1.884e17 s - 6.722e19
   5:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                                -4.212e09 s - 3.132e14
   6:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
Continuous-time transfer function.

step(G, 0.02)

dcgain(G)

ans = 6×2

1.0000 0.0002 -0.0000 0 Inf Inf 1.0000 -0.0000 -0.0000 -9.4309 0.0028 -0.0000

Gtest = G(3,1)      % transfer function of Xi(s)/U(s)

Gtest =
 
     s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.044e13 s^2 + 1.998e16 s + 537.4
  ----------------------------------------------------------------------------
  s^6 + 1.062e05 s^5 + 3.194e09 s^4 + 3.045e13 s^3 + 2.065e16 s^2 + 7.127e18 s
 
Continuous-time transfer function.

[num, den] = tfdata(Gtest, 'v')

num = 1×7

1.0e+16 * 0 0.0000 0.0000 0.0000 0.0030 1.9980 0.0000

den = 1×7

1.0e+18 * 0.0000 0.0000 0.0000 0.0000 0.0207 7.1274 0

num(7)

ans = 537.4219

Actually, I got a different value on an older machine, but I still same output as Paul did, at least up to 0.02 sec.

[John]

@Sam Chak Thanks. Interesting: I do get the same TFs that you showed. My numeric outputs match yours, specifically the 3rd state from input 1 to output (, u = ref) are the same. My plots also match when using G.

Except, when I pass the SS system into step() or impulse() -- then there are 0s for the 3rd row.

Check out the attached code that compares G and sysCl_lqgi (uses the same .mat file). It's also pasted below to Run in this post.

Here's what I see:

1) On my local machine, see the plot immediately below: xi = 0s for the SS sys, and nonzero for G. That seems odd, since G was made from the SS sys.

2) But when the code is Run locally in this post, all plots non-0 and match as expected (very bottom of post)

Also, Paul and I both ran the same code with the same variable, with the same output on Answers console.

But, when I run it locally on my machine after clearing everything, i do not (seem to) get the same output, and only for that variable.

Are you able to locally run the code (actually your code :) ) I added here, with the .mat file variable? If you run that code, what do you get for the plot?

clear all
load('sys.mat')
A   = sysCl_lqgi.A;
B   = sysCl_lqgi.B;
C   = eye(6);
D   = 0;
sys = ss(A, B, C, D);
G   = tf(sys)

G =
 
  From input 1 to output...
                         1.596e10 s^2 + 6.747e14 s + 7.127e18
   1:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                 1.596e10 s^3 + 6.747e14 s^2 + 7.127e18 s - 5.948e-11
   2:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
          s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.044e13 s^2 + 1.998e16 s + 537.4
   3:  ----------------------------------------------------------------------------
       s^6 + 1.062e05 s^5 + 3.194e09 s^4 + 3.045e13 s^3 + 2.065e16 s^2 + 7.127e18 s
 
                         1.596e10 s^2 + 6.747e14 s + 7.127e18
   4:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                 1.596e10 s^3 + 6.747e14 s^2 + 7.127e18 s - 5.564e06
   5:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
              s^4 + 1.062e05 s^3 + 3.194e09 s^2 + 3.044e13 s + 1.998e16
   6:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
  From input 2 to output...
                    1e05 s^3 + 1.06e10 s^2 + 3.177e14 s + 1.596e15
   1:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                  1e05 s^4 + 1.06e10 s^3 + 3.177e14 s^2 + 1.596e15 s
   2:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                     -1e05 s^3 - 1.06e10 s^2 - 3.177e14 s - 1.596e15
   3:  ----------------------------------------------------------------------------
       s^6 + 1.062e05 s^5 + 3.194e09 s^4 + 3.045e13 s^3 + 2.065e16 s^2 + 7.127e18 s
 
                          4.212e09 s^2 + 3.132e14 s - 0.2129
   4:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                         4.398e13 s^2 - 1.884e17 s - 6.722e19
   5:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                                -4.212e09 s - 3.132e14
   6:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
Continuous-time transfer function.

[y1, t1, x1] = step(G, 0.02);
[y2, t2, x2] = step(sysCl_lqgi, 0.02);

dcgain(G)

ans = 6×2

1.0000 0.0002 -0.0000 0 Inf Inf 1.0000 -0.0000 -0.0000 -9.4309 0.0028 -0.0000

dcgain(sysCl_lqgi)

ans = 1×2

1.0000 0.0002

Gtest = G(3,1) 

Gtest =
 
     s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.044e13 s^2 + 1.998e16 s + 537.4
  ----------------------------------------------------------------------------
  s^6 + 1.062e05 s^5 + 3.194e09 s^4 + 3.045e13 s^3 + 2.065e16 s^2 + 7.127e18 s
 
Continuous-time transfer function.

[num, den] = tfdata(Gtest, 'v')

num = 1×7

1.0e+16 * 0 0.0000 0.0000 0.0000 0.0030 1.9980 0.0000

den = 1×7

1.0e+18 * 0.0000 0.0000 0.0000 0.0000 0.0207 7.1274 0

num(7)

ans = 537.4219

figure
subplot(2, 1, 1);
plot(t1, y1(:, 1, 1), 'b', 'LineWidth', 1.5);
hold on
plot(t2, x2(:, 1, 1), 'k:', 'LineWidth', 2.5)
grid on;
title('x1', 'FontSize', 14)
legend({'$x_1$ G', '$x_1$ Sys'},'Interpreter','latex', 'FontSize', 16);
xlim([0, 0.02]);

subplot(2, 1, 2);
plot(t1, y1(:, 3, 1), 'b', 'LineWidth', 1.5);
hold on
plot(t2, x2(:, 3, 1), 'k:', 'LineWidth', 2.5);
grid on;
title('xi', 'FontSize', 14)
legend({'$x_i G$', '$x_i$ Sys'},'Interpreter','latex', 'FontSize', 16);
xlim([0, 0.02]);
xlabel('Time [s]', 'FontSize', 14)
sgtitle('STEP(): Ref response');

[Walter Roberson]

/Applications/MATLAB_R2022b.app/toolbox/shared/controllib/engine/@DynamicSystem/step.m % ss method

is the R2022b output on a MacOS system for the same source as

/MATLAB/toolbox/shared/controllib/engine/@DynamicSystem/step.m % ss method

The only difference is the installation directory prefix, /MATLAB in one case an /Applications/MATLAB_R2022b.app in the other case.

[John]

@Sam Chak @Paul Check out this modification to the first part of the code, with the option switch.

Option 3 works correctly on my machine, and option 1 and 2 do not.

Do all 3 plot correctly for you?

Even if there's something wrong with my SS espression, it seems unexpected that @Paul and my outputs would be different with the same code and same variable.

clear all
load('sys.mat')
A   = sysCl_lqgi.A;
B   = sysCl_lqgi.B;
C   = eye(6);
D   = 0;
sys = ss(A, B, C, D);
G   = tf(sys)

G =
 
  From input 1 to output...
                         1.596e10 s^2 + 6.747e14 s + 7.127e18
   1:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                 1.596e10 s^3 + 6.747e14 s^2 + 7.127e18 s - 5.948e-11
   2:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
          s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.044e13 s^2 + 1.998e16 s + 537.4
   3:  ----------------------------------------------------------------------------
       s^6 + 1.062e05 s^5 + 3.194e09 s^4 + 3.045e13 s^3 + 2.065e16 s^2 + 7.127e18 s
 
                         1.596e10 s^2 + 6.747e14 s + 7.127e18
   4:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                 1.596e10 s^3 + 6.747e14 s^2 + 7.127e18 s - 5.564e06
   5:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
              s^4 + 1.062e05 s^3 + 3.194e09 s^2 + 3.044e13 s + 1.998e16
   6:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
  From input 2 to output...
                    1e05 s^3 + 1.06e10 s^2 + 3.177e14 s + 1.596e15
   1:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                  1e05 s^4 + 1.06e10 s^3 + 3.177e14 s^2 + 1.596e15 s
   2:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                     -1e05 s^3 - 1.06e10 s^2 - 3.177e14 s - 1.596e15
   3:  ----------------------------------------------------------------------------
       s^6 + 1.062e05 s^5 + 3.194e09 s^4 + 3.045e13 s^3 + 2.065e16 s^2 + 7.127e18 s
 
                          4.212e09 s^2 + 3.132e14 s - 0.2129
   4:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                         4.398e13 s^2 - 1.884e17 s - 6.722e19
   5:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
                                -4.212e09 s - 3.132e14
   6:  ------------------------------------------------------------------------
       s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.045e13 s^2 + 2.065e16 s + 7.127e18
 
Continuous-time transfer function.

option = 3;
switch option
    case 1
        % Bad: xi = 0
        C = [1, zeros(1, 5)];
        D = 0;
    case 2
        % Bad: xi = 0
        C = zeros(6);
        C(1) = 1;
        D = zeros(6,2);
    case 3
        % Good: xi matches
        C = eye(6);
        D = zeros(6,2);
    otherwise
        disp('error');
        return
end
% Make different versions of sysCl_lqgi
sysCl_lqgi = ss(sysCl_lqgi.A, sysCl_lqgi.B, C, D)

sysCl_lqgi =
 
  A = 
               x1          x2          x3          x4          x5          x6
   x1           0           1           0           0           0           0
   x2  -6.317e+04        -160           0  -4.468e+07  -6.376e+04   1.596e+10
   x3          -1           0           0           0           0           0
   x4   4.212e+04           0           0  -4.212e+04           1           0
   x5   4.398e+08           0           0  -4.845e+08  -6.392e+04   1.596e+10
   x6           0           0           0          -1           0           0
 
  B = 
          u1     u2
   x1      0      0
   x2      0  1e+05
   x3      1      0
   x4      0      0
   x5      0      0
   x6      1      0
 
  C = 
       x1  x2  x3  x4  x5  x6
   y1   1   0   0   0   0   0
   y2   0   1   0   0   0   0
   y3   0   0   1   0   0   0
   y4   0   0   0   1   0   0
   y5   0   0   0   0   1   0
   y6   0   0   0   0   0   1
 
  D = 
       u1  u2
   y1   0   0
   y2   0   0
   y3   0   0
   y4   0   0
   y5   0   0
   y6   0   0
 
Continuous-time state-space model.

[y1, t1, x1] = step(G, 0.02);
[y2, t2, x2] = step(sysCl_lqgi, 0.02);

dcgain(G)

ans = 6×2

1.0000 0.0002 -0.0000 0 Inf Inf 1.0000 -0.0000 -0.0000 -9.4309 0.0028 -0.0000

dcgain(sysCl_lqgi)

ans = 6×2

1.0000 0.0002 0 0 Inf Inf 1.0000 0 -0.0000 -9.4309 0.0028 -0.0000

Gtest = G(3,1) 

Gtest =
 
     s^5 + 1.062e05 s^4 + 3.194e09 s^3 + 3.044e13 s^2 + 1.998e16 s + 537.4
  ----------------------------------------------------------------------------
  s^6 + 1.062e05 s^5 + 3.194e09 s^4 + 3.045e13 s^3 + 2.065e16 s^2 + 7.127e18 s
 
Continuous-time transfer function.

[num, den] = tfdata(Gtest, 'v')

num = 1×7

1.0e+16 * 0 0.0000 0.0000 0.0000 0.0030 1.9980 0.0000

den = 1×7

1.0e+18 * 0.0000 0.0000 0.0000 0.0000 0.0207 7.1274 0

num(7)

ans = 537.4219

figure
subplot(2, 1, 1);
plot(t1, y1(:, 1, 1), 'b', 'LineWidth', 1.5);
hold on
plot(t2, x2(:, 1, 1), 'k:', 'LineWidth', 2.5)
grid on;
title('x1', 'FontSize', 14)
legend({'$x_1$ G', '$x_1$ Sys'},'Interpreter','latex', 'FontSize', 16);
xlim([0, 0.02]);

subplot(2, 1, 2);
plot(t1, y1(:, 3, 1), 'b', 'LineWidth', 1.5);
hold on
plot(t2, x2(:, 3, 1), 'k:', 'LineWidth', 2.5);
grid on;
title('xi', 'FontSize', 14)
legend({'$x_i G$', '$x_i$ Sys'},'Interpreter','latex', 'FontSize', 16);
xlim([0, 0.02]);
xlabel('Time [s]', 'FontSize', 14)
sgtitle('STEP(): Ref response');

[John]

@Walter Roberson Thanks for verifying. I thought it was the same directory as well, ie it seems like the matching function is called.

[John]

@Paul @Sam Chak @Walter Roberson Just to update as requested, this is the response from Mathworks so far:

"Thank you for your patience in awaiting my reply. I was able to reproduce this issue. I am unsure as to the cause of this. The error arises in R2022b and all previous versions but disappears in R2023a onwards. The only ostensible change in "step" in R2023a is that the default start time is not limited to t = 0. This does not apply to your script, since the time array started at 0.

Some of my colleagues and I are taking a look at the moment. I will get back to you when I have further information to share."

[Sam Chak]

Hi @John, thanks for your update. It's good that you pointed out the issue. Luckily in my work, I seldom use high-value parameters.