LQI controller is weirdly slow in response and doesn't always reach a goal on plant

EDIT: Here is data I used for identification

I am trying to develop a servomechanism LQI controller. The servo is supposed to act as an actuator for a brake pedal - it is connected to a pedal via steel wire and therefore pushes on the pedal similarily to how a driver would. My goal was to first implement the control in MATLAB step by step and then rewrite it to embedded C. In simulation the controller behaves in a very satisfactory way, but when moving to plant, the response seems very delayed. Additionaly, in simulation I have neved had any issues with goal reaching, but on plant it turned out that the maximum the controller can achieve is 10 bars, and even when the goal pressure is higher, it never reaches it, just settles on the lower value. What is weird is that there is no issue in reaching lower pressures. I have tried to change the integral gain value since in the very beginning of testing, the maximum reachable pressure was even lower, and somehow changing the gain got me up to 10 bars - but never more. I am new to control theory. I need this system to be way faster and I want to make the delay the smallest possible, preferably around 30 ms.

Below I'm adding two pictures - one is from simulation and the other one directly from the plant. The delay on the plant is around 1 second (around 1100 ms).

% System identification
A = [0.748153090242338 -0.0770817542371630 -0.272770107988345 -0.0680287168856025 0.274040229219026;
    0.0264773037738115 0.802405697397900 -0.244152268701560 0.0743061296316669 -0.0764153326240192;
    -0.0409228451430832 0.360565865830120 0.715237466131514 0.0837023187613815 0.199733910116335;
    0.0884449413535669 -0.0520801836735802 -0.0862407832241391 0.982500599129848 -0.0782458817344473;
    -0.00176167556838596 0.145291619650476 -0.179472954855322 -0.156570713297809 0.410924513958431];

B = [-0.0691336950664988;
    0.298054008136356;
    -0.101282668881553;
    -0.135571075515549;
    -0.197136001368224];

C = [0.395474090810162 0.0874731850181924 0.127527899108947 -0.0511706805106379 -0.0420611557807524];

D = 0;

sys = ss(A, B, C, D, 0.02);
step(sys);

[A_size_row, A_size_column] = size(A);
Q = 70*eye(A_size_row+1);
R = 5;

A_aug = [A, zeros(A_size_column, 1); -C, 0];
B_aug = [B; 0];

[K, ~, ~] = lqr(A_aug, B_aug, Q, R); % LQR solver
K_x = K(1:5);
K_i = K(6);

dt = 0.02;  % time step;
T = 1000;  % sim time in milliseconds (or is it seconds?)
time = 0:dt:T;

x = [0; 
    0;
    0;
    0;
    0];  % initial state
x_i = 0;  % initial integrator state
y = [];  % output storage
u = [];  % control effort storage
ref = [];  % reference signal storage
states = []; % state storage
r=0;

% Simulation loop
for t = time
    % Reference signal
    if t == 50
        r = 5;
    end
    
    % Control law
    u_t = -K_x * x - K_i * x_i;

    % System dynamics
    x_dot = A * x + B * u_t;
    x = x + x_dot * dt;
    
    % Integrator dynamics
    y_t = C * x;
    x_i_dot = r - y_t;
    x_i = x_i + x_i_dot * dt;
    
    % Store results
    y = [y, y_t];
    u = [u, u_t];
    ref = [ref, r];
    states = [states, x]; % Store state vector
end

% Plot results
figure;
subplot(3,1,1);
plot(time, ref, 'r--', 'LineWidth', 1.5);
hold on;
plot(time, y, 'b', 'LineWidth', 1.5);
xlabel('Time (ms)');
ylabel('Output y');
title('System Output vs Reference Signal');
legend('Reference', 'Output');
grid on;

subplot(3,1,2);
plot(time, u, 'g', 'LineWidth', 1.5);
xlabel('Time (ms)');
ylabel('Control Effort u');
title('Control Effort');
grid on;

subplot(3,1,3);
plot(time, states, 'LineWidth', 1.5);
xlabel('Time (ms)');
ylabel('States');
title('States over Time');
legend('x1', 'x2', 'x3', 'x4', 'x5');
grid on;

[Mathieu NOE]

hello

maybe the simulation gives you good results because the LQI regulator is well designed vs the model you give him - this trigger my question , how sure are you about your model ? - looks like a first oder + delay , which BTW could probably be controlled with some simpler PID controller ?

also the 200 ms delay in your ss model looks quite large and again this makes me believe there is probably a issue between reality and model. How can your actuator have such dead time ??

second as your LQI control law is more difficult to implement (vs PID) in real time , maybe there is a bug somewhere ... did you compare the C code vs the matlab code before real time testing ?

[Wiktoria]

Hello!

Thank you for your answer. There might be some issue about the model - I used system identification toolbox and tried to follow MATLAB guides about how to do that. After I got the response to the step from the real time plant, I tried to fit it with the toolbox using the "estimate state-space models" option.

The actuator itself doesn't have such a big delay at all! Actually, the delay is around 20 ms only. When I set the pulse directly ("by hand") it's very fast, but when controller does that - it suddenly is very slow. So you're probably right, it might be the issue with the model. I compared the code from matlab before testing it on plant - it seemed identical.

My goal is to use LQI or LQR as I heard it's more precise and robust, as well as I want to simply try to implement it to gain knowledge (and it's also my project that I need to finish). Is there any way to make the response faster or to fit the proper, faster model? Thank you for your help!

[Paul]

Hi Wiktoria,

These lines define the system with a discrete-time model:

A = [0.748153090242338 -0.0770817542371630 -0.272770107988345 -0.0680287168856025 0.274040229219026;
    0.0264773037738115 0.802405697397900 -0.244152268701560 0.0743061296316669 -0.0764153326240192;
    -0.0409228451430832 0.360565865830120 0.715237466131514 0.0837023187613815 0.199733910116335;
    0.0884449413535669 -0.0520801836735802 -0.0862407832241391 0.982500599129848 -0.0782458817344473;
    -0.00176167556838596 0.145291619650476 -0.179472954855322 -0.156570713297809 0.410924513958431];
B = [-0.0691336950664988;
    0.298054008136356;
    -0.101282668881553;
    -0.135571075515549;
    -0.197136001368224];
C = [0.395474090810162 0.0874731850181924 0.127527899108947 -0.0511706805106379 -0.0420611557807524];
D = 0;
sys = ss(A, B, C, D, 0.02);
isdt(sys)

ans = logical
   1

But these lines

A_aug = [A, zeros(A_size_column, 1); -C, 0];
B_aug = [B; 0];

are what would be used to augment a continuous-time plant with an integrator.

And this form of the call to lqr

[K, ~, ~] = lqr(A_aug, B_aug, Q, R); % LQR solver

only applied for continuous-time systems. So it appears there is a disconnect between the form of the model and how the control law is being designed.

If you want to use lqr, you'll need to augment the plant with a discrete-time integrator and call lqr with an ss object as the first argument. Alternatively, you can use lqi to design the gains.

Finally, the simulation loop is treating the state space model as if it describes the system dynamics in continuous time, even though sys is a discrete-time model.

Before doing anything else the actual domain (continuous or discrete) of A,B,C,D needs to be established.

[Wiktoria]

Hello, thank you for your answer!

I have checked that the augmentation for discrete-time plant is supposedly

where m is a number of states. I don't quite understand if it is 5x5 matrix? But from what I understand, if A is 5x5 matrix, -C is 1x5, therefore zeros must be 5x1, so I matrix should be 1x1?

A_aug = [A, zeros(A_size_column, 1); -C, 1];
B_aug = [B; 0];

Aufter augmenting I computed the gains as you mentioned, using the formula below:

sys = ss(A, B, C, D, 0.02);
[K2,S,e] = lqi(sys,Q,R);

After generating a step with a closed loop system matrix, it turned out it looked pretty good.

Acl=A-B*K(1:5);
sys_lqr=ss(Acl,B,C,D, 0.02);
step(sys_lqr);

I am confused though, you mentioned that the simulation loop is treating the state space model as if it describes the system dynamics in continuous time, even though sys is a discrete-time model. Could you elaborate what you mean? Or is there perhaps any way to get to know how it should look like for discrete time? I tried to use the gains from lqi in the loop I have now (so the one for continuous time) and it is unstable (which isn't really surprising). How should the loop look like? I am new to control theory and didn't know there is a difference in how the loop is treating the space model - is there a difference in equations? How to correct it?

[Paul]

Let's go back to step 1. You have a model with A,B,C,D matrices defined as

A = [0.748153090242338 -0.0770817542371630 -0.272770107988345 -0.0680287168856025 0.274040229219026;
    0.0264773037738115 0.802405697397900 -0.244152268701560 0.0743061296316669 -0.0764153326240192;
    -0.0409228451430832 0.360565865830120 0.715237466131514 0.0837023187613815 0.199733910116335;
    0.0884449413535669 -0.0520801836735802 -0.0862407832241391 0.982500599129848 -0.0782458817344473;
    -0.00176167556838596 0.145291619650476 -0.179472954855322 -0.156570713297809 0.410924513958431];
B = [-0.0691336950664988;
    0.298054008136356;
    -0.101282668881553;
    -0.135571075515549;
    -0.197136001368224];
C = [0.395474090810162 0.0874731850181924 0.127527899108947 -0.0511706805106379 -0.0420611557807524];
D = 0;

Do these matrices represent a model of the form:

dx/dt = A*x(t) + B*u(t)

y(t) = C*x(t) + D*u(t)

or of the form:

x(k+1) = A*x(k) + B*u(k)

y(k) = C*x(k) + D*u(k)

?

If the former, what was the intent of the 0.02 parameter in the line

sys = ss(A, B, C, D, 0.02);

If the latter, is the 0.02 supposed to be the sampling period? Like maybe you developed the A,B,C,D matrices using system indentification on data that was sampled every 0.02 seconds?

[Wiktoria]

It was the latter, I wanted to represent a form

The 0.02s is supposed to be the sampling period, yes. I have developed the system using System Identification Toolbox using data that was sampled every 0.02s. In the toolbox, I was supposed to input the sampling period, so that's what I did.

[Paul]

Discussion continued at this answer.

[Pavl M.]

Attn.:

The original idea to employ the Hamiltonian matrix for the specific closed loop agent design and simulation:

Acl=A-B*K; % which is H matrix in my answer.

sys_lqr=ss(Acl,B,C,D, s_time);

step(sys_lqr);

is my, not of that Paul, please let's keep order who is who.

It is really good question.

Kindly see my recent the solution update:

(Updated, 17 Nov. 2024, see chk1.m script, Final for the numerical LQR, LQI gains design):

Here is a solution for much more faster LQ controller(agent) working in loop ( see my answer code + m going to add simulation of your loop with noise addition) with reduced oscillations by pre-filter I developed(constructed) novel (plagiarism free, not taken from any published manuscript, whatever indeed works in soft real time embedded). What is better for you LQR(no augmentation of matrices A,B,C,D, is required) or LQI(needs augmentation of A,B,C,D, matrices)Carefull selection of Q,R,N matrices are similar in goal vs the process of PID tuning.

Let's see maybe you can use the controller (smart agent for which the good K gains finding we have completed) with your system without complex Particle Filter or Kalman state X estimator if to rely on your linearized matrices A,B,C,D (please check are you sure about the matrices, and what is the actual physical system to control in non-linear, time varying case? whether they work in embedded code?) solely then the state X can be computed inside the embedded code from equation

The control u computed on your embedded board relying on the actual real physical system estimation matrices A,B,C,D, hard coded in your embedded code and using actual physical system output y, the control u computed you just connect via appropriate transducer (signal conditioner) to your actual physical system actuator (servo motor or any other valid actuator).

Hope to construct together many or a few most valuable, really utile, useful actor controllers(agents) embeddable low code and high code C/C++ embedded running on digital computer codes for industry. Let's focus on this workflow.

Answer 1

[Mathieu NOE]

hello again

used a simplified version of your model and applied simple PID control (easier for real time conversion) :

% Define system parameters
K = 1;  % Gain
T = 0.1;  % Time constant
L = 0.2;  % Delay

% Transfer function of the first-order system with delay
s = tf('s');
sys_delayed = exp(-L*s)*tf(K, [T, 1]);



% Define PID controller gains
Kp = 0.6;
Ki = 5*Kp;
Kd = 0;

% Transfer function of the PID controller
pid = tf([Kd, Kp, Ki], [1, 0]);

% Closed-loop system
closed_loop_sys = feedback(pid*sys_delayed,1);

% Step response
t = linspace(0, 5, 1000);
[yOL, t] = step(sys_delayed, t);
[yCL, t] = step(closed_loop_sys, t);

% Plot the step response
plot(t, yOL,'b',t, yCL,'r')
xlabel('Time (s)')
ylabel('Response')
title('Step Response of First-Order System with Delay and PID Controller')
grid on

[Mathieu NOE]

a further attemp to optimize the PI controller using pidtune

% Define system parameters
K = 1;  % Gain
T = 0.1;  % Time constant
L = 0.2;  % Delay
% Transfer function of the first-order system with delay
s = tf('s');
sys_delayed = exp(-L*s)*tf(K, [T, 1]);

% Define PID controller gains
[pid,info] = pidtune(sys_delayed,'PI');

% Closed-loop system
closed_loop_sys = feedback(pid*sys_delayed,1);
% Step response
t = linspace(0, 5, 1000);
[yOL, t] = step(sys_delayed, t);
[yCL, t] = step(closed_loop_sys, t);
% Plot the step response
plot(t, yOL,'b',t, yCL,'r')
xlabel('Time (s)')
ylabel('Response')
title('Step Response of First-Order System with Delay and PID Controller')
grid on

[Pavl M.]

that is too simple, trivial, simplification often makes losses in features and original beauty of complexity, make it simple, but not simpler. In this question case the problem is another, it is of proper LQ controller design (the wise cost matrices Q,R,N selection and state estimation or better recovery LTM) and LQ-derivatives compare-contrast and most fitted selection implementable on embedded board must be original the 5+ size A,B,C,D, matrices and not what you/others pushed, also the delay and amplitude defficiency they did not tackle, while I handle in my answers/solution pathes, LQ may be better in performance than PID since it inherently deals with convex optimization and I found the gain K matrix for reduced overshoot, faster settling, more close to desired gain and lower rise time (and so less sluggish) as per the original requirements. Hire me for better, including block diagrams/block schemes and flow of control with TCE Simulink I can make normal.

[Mathieu NOE]

well the art of engineering is to choose the right tool for the right job. I have nothing against the beauty of complexity, but afetr many yeras in the industry, 99.9% of the time I prefer to simplify the problem first before designing a solution for it. of course , it's everyone personal choice / preferences to pick this or that appoach.

but when it comes down to embedded software, the complexier your design the more headaches to test and make sure sure there is no hidden bug somewhere.

Answer 2

[Pavl M.]

It is indeed really factually very valued, yet undervalued and quite intersting, utile question.

(Updated, 21 Nov. 2024, Final for the numerical LQR, LQI gains design):

Here is a solution for much more faster LQ controller(agent) working in loop with reduced oscillations by pre-filter I developed(constructed) novel (plagiarism free, not taken from any published manuscript, further if will be found interested engineering specialist to hire me I can analyze and make loop transfer recovery (LTR) and LTM agents(regulators,controllers) and tackle other sophisticated and simple looped systems, whatever indeed works in soft real time embedded). What is better for you LQR(no augmentation of matrices A,B,C,D, is required) or LQI(needs augmentation of A,B,C,D, matrices), see how TCE NCE MPPL Matlab internal lqr state feedback controlle works faster without augmentation, to where the setpoint comes to just state feedback built-in lqr:

clc
clear all
close all

s_time = 0.02; %[sec], SI by default in TCE Matlab
Np = 0.01; %noise power
nstates = 5;
ninputs = 1;
ncontrols = 1;
noutputs = 1;

% System identification
A = [0.748153090242338 -0.0770817542371630 -0.272770107988345 -0.0680287168856025 0.274040229219026;
    0.0264773037738115 0.802405697397900 -0.244152268701560 0.0743061296316669 -0.0764153326240192;
    -0.0409228451430832 0.360565865830120 0.715237466131514 0.0837023187613815 0.199733910116335;
    0.0884449413535669 -0.0520801836735802 -0.0862407832241391 0.982500599129848 -0.0782458817344473;
    -0.00176167556838596 0.145291619650476 -0.179472954855322 -0.156570713297809 0.410924513958431];
B = [-0.0691336950664988;
    0.298054008136356;
    -0.101282668881553;
    -0.135571075515549;
    -0.197136001368224];
C = [0.395474090810162 0.0874731850181924 0.127527899108947 -0.0511706805106379 -0.0420611557807524];
D = 0;

sys1 = ss(A, B, C, D, s_time);
sys2 = ss(A,B,C,D);

disp('Whether the system is stable, minimum phase and proper:')

Whether the system is stable, minimum phase and proper:

isstable(sys1)

ans = logical
   1

isminphase(tfdata(tf(sys1),'v'))

ans = logical
   0

isproper(sys1)

ans = logical
   1

Qc= ctrb(A,B)

Qc = 5×5

-0.0691 -0.0919 -0.0764 -0.0676 -0.0674 0.2981 0.2670 0.2051 0.1271 0.0493 -0.1013 -0.0129 0.0803 0.1345 0.1429 -0.1356 -0.1307 -0.1495 -0.1761 -0.2023 -0.1971 0.0018 0.0625 0.0646 0.0486

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

controllab = rank(Qc)

controllab = 5

x = 1;
disp('The pair (A,B) must be stabilizable')

The pair (A,B) must be stabilizable

ic = [0 0 0 0 0];

figure
step(sys1);

stepinfo(sys1)

ans = struct with fields:
         RiseTime: 0.2200
    TransientTime: 0.5600
     SettlingTime: 0.5600
      SettlingMin: 0.7101
      SettlingMax: 0.7932
        Overshoot: 1.8688
       Undershoot: 1.3344
             Peak: 0.7932
         PeakTime: 0.7000

hasdelay(sys1)

ans = logical
   0

td = totaldelay(sys1)

td = 0

sysnd = absorbDelay(sys1)

sysnd =
 
  A = 
              x1         x2         x3         x4         x5
   x1     0.7482   -0.07708    -0.2728   -0.06803      0.274
   x2    0.02648     0.8024    -0.2442    0.07431   -0.07642
   x3   -0.04092     0.3606     0.7152     0.0837     0.1997
   x4    0.08844   -0.05208   -0.08624     0.9825   -0.07825
   x5  -0.001762     0.1453    -0.1795    -0.1566     0.4109
 
  B = 
             u1
   x1  -0.06913
   x2    0.2981
   x3   -0.1013
   x4   -0.1356
   x5   -0.1971
 
  C = 
             x1        x2        x3        x4        x5
   y1    0.3955   0.08747    0.1275  -0.05117  -0.04206
 
  D = 
       u1
   y1   0
 
Sample time: 0.02 seconds
Discrete-time state-space model.

sysx = pade(sys1,nstates)

Warning: The model is already discrete. For discrete models, the "pade" command maps all delays to poles at z=0.

sysx =
 
  A = 
              x1         x2         x3         x4         x5
   x1     0.7482   -0.07708    -0.2728   -0.06803      0.274
   x2    0.02648     0.8024    -0.2442    0.07431   -0.07642
   x3   -0.04092     0.3606     0.7152     0.0837     0.1997
   x4    0.08844   -0.05208   -0.08624     0.9825   -0.07825
   x5  -0.001762     0.1453    -0.1795    -0.1566     0.4109
 
  B = 
             u1
   x1  -0.06913
   x2    0.2981
   x3   -0.1013
   x4   -0.1356
   x5   -0.1971
 
  C = 
             x1        x2        x3        x4        x5
   y1    0.3955   0.08747    0.1275  -0.05117  -0.04206
 
  D = 
       u1
   y1   0
 
Sample time: 0.02 seconds
Discrete-time state-space model.

[DelayFreeSys,tau] = getDelayModel(sys1)

DelayFreeSys =
 
  A = 
              x1         x2         x3         x4         x5
   x1     0.7482   -0.07708    -0.2728   -0.06803      0.274
   x2    0.02648     0.8024    -0.2442    0.07431   -0.07642
   x3   -0.04092     0.3606     0.7152     0.0837     0.1997
   x4    0.08844   -0.05208   -0.08624     0.9825   -0.07825
   x5  -0.001762     0.1453    -0.1795    -0.1566     0.4109
 
  B = 
             u1
   x1  -0.06913
   x2    0.2981
   x3   -0.1013
   x4   -0.1356
   x5   -0.1971
 
  C = 
             x1        x2        x3        x4        x5
   y1    0.3955   0.08747    0.1275  -0.05117  -0.04206
 
  D = 
       u1
   y1   0
 
Sample time: 0.02 seconds
Discrete-time state-space model.

tau =

  0x1 empty double column vector

Gd1 = c2d(sys2,s_time,'impulse')

Gd1 =
 
  A = 
               x1          x2          x3          x4          x5
   x1       1.015   -0.001577   -0.005541   -0.001399    0.005536
   x2   0.0005412       1.016   -0.004958    0.001511   -0.001557
   x3  -0.0008271    0.007327       1.014    0.001702    0.004031
   x4      0.0018    -0.00107   -0.001754        1.02   -0.001585
   x5  -3.619e-05     0.00293   -0.003634   -0.003176       1.008
 
  B = 
              u1
   x1   -0.00142
   x2   0.006069
   x3   -0.00203
   x4  -0.002764
   x5  -0.003942
 
  C = 
             x1        x2        x3        x4        x5
   y1    0.3955   0.08747    0.1275  -0.05117  -0.04206
 
  D = 
              u1
   y1  2.088e-05
 
Sample time: 0.02 seconds
Discrete-time state-space model.

Gd2 = d2d(sys1,s_time,'tustin')

Gd2 =
 
  A = 
              x1         x2         x3         x4         x5
   x1     0.7482   -0.07708    -0.2728   -0.06803      0.274
   x2    0.02648     0.8024    -0.2442    0.07431   -0.07642
   x3   -0.04092     0.3606     0.7152     0.0837     0.1997
   x4    0.08844   -0.05208   -0.08624     0.9825   -0.07825
   x5  -0.001762     0.1453    -0.1795    -0.1566     0.4109
 
  B = 
             u1
   x1  -0.06913
   x2    0.2981
   x3   -0.1013
   x4   -0.1356
   x5   -0.1971
 
  C = 
             x1        x2        x3        x4        x5
   y1    0.3955   0.08747    0.1275  -0.05117  -0.04206
 
  D = 
       u1
   y1   0
 
Sample time: 0.02 seconds
Discrete-time state-space model.

%figure
%step(Gd1)
figure
step(Gd2)

[A_size_row, A_size_column] = size(A);
Q = 70*eye(A_size_row+1)

Q = 6×6

70 0 0 0 0 0 0 70 0 0 0 0 0 0 70 0 0 0 0 0 0 70 0 0 0 0 0 0 70 0 0 0 0 0 0 70

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

R = 7 %*eye(ninputs);

R = 7

R2 = 5 %*eye(ninputs+1);

R2 = 5

N = ones(nstates,1)

N = 5×1

1 1 1 1 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

N2 = ones(nstates+1,1)

N2 = 6×1

1 1 1 1 1 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%S
%E
%P %matrices to find furthermore.
A_aug = [A, zeros(A_size_column, 1); -C, 0];
B_aug = [B; 0];
Q2 = 30*eye(A_size_row) % To be selected using Adaptive Dynamic Programming or 

Q2 = 5×5

30 0 0 0 0 0 30 0 0 0 0 0 30 0 0 0 0 0 30 0 0 0 0 0 30

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%similar heuristic for rubustness
%R = 5;
%S = eye(nstates,noutputs);
%E = eye(nstates,nstates);
%[~,p] = chol(R)
%% Just ARE
[X1, L1, G1] = care(A_aug, B_aug, Q, [1]) %, S, E)

X1 = 6×6

1.0e+07 * 0.3282 -0.1878 -0.3501 -0.5924 0.1873 -0.0026 -0.1878 0.1265 0.2121 0.3805 -0.1127 0.0010 -0.3501 0.2121 0.3817 0.6571 -0.2033 0.0024 -0.5924 0.3805 0.6571 1.1607 -0.3502 0.0037 0.1873 -0.1127 -0.2033 -0.3502 0.1086 -0.0013 -0.0026 0.0010 0.0024 0.0037 -0.0013 0.0001

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

L1 = 
  -0.1512 + 0.0000i
  -0.6857 + 0.0000i
  -0.7866 + 0.2586i
  -0.7866 - 0.2586i
  -0.8789 + 0.0000i
  -3.3905 + 0.0000i

G1 = 1×6

1.0e+03 * 2.0157 -1.7673 -2.4339 -5.1361 1.3512 -0.0084

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%% For Just LQR:
[K1,S1,P1] = lqr(sys1,Q2,R)

K1 = 1×5

-0.5461 1.3958 0.1922 -1.2207 -0.1025

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

S1 = 5×5

77.7954 -18.7303 -48.1839 34.8112 2.5513 -18.7303 111.1492 33.6975 11.7969 10.2405 -48.1839 33.6975 114.2713 -47.1588 7.3914 34.8112 11.7969 -47.1588 321.9632 -53.5409 2.5513 10.2405 7.3914 -53.5409 50.5475

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

P1 = 
   0.6870 + 0.2420i
   0.6870 - 0.2420i
   0.6773 + 0.0000i
   0.6323 + 0.0000i
   0.3556 + 0.0000i

%% For just LQI:
[Klqi,Slqi,elqi] = lqi(sys1,Q,R2,N2)

Klqi = 1×6

-0.6764 1.8167 0.3253 -1.7091 -0.0709 -1.5924

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Slqi = 6×6

1.0e+04 * 0.0185 -0.0034 -0.0111 0.0054 0.0011 -0.0356 -0.0034 0.0225 0.0073 0.0056 0.0024 0.0050 -0.0111 0.0073 0.0265 -0.0101 0.0016 0.0059 0.0054 0.0056 -0.0101 0.0782 -0.0139 0.0953 0.0011 0.0024 0.0016 -0.0139 0.0122 -0.0258 -0.0356 0.0050 0.0059 0.0953 -0.0258 1.8142

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

elqi = 
   0.9959 + 0.0000i
   0.6731 + 0.2256i
   0.6731 - 0.2256i
   0.6576 + 0.0531i
   0.6576 - 0.0531i
   0.2009 + 0.0000i

%In discrete time, lqi computes the integrator output xi 
%using the forward Euler formula
%xi[n+1]=xi[n]+Ts(r[n]−y[n])
%where Ts is the sample time of SYS.
%What if to try trapezoidal and Backward-Euler integration methods?
[K2,S2,P2] = dlqr(A,B,Q2,R,N)

K2 = 1×5

-0.5260 1.3925 0.2548 -1.2030 -0.0769

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

S2 = 5×5

79.0469 -19.5596 -47.6365 36.6696 3.5117 -19.5596 106.8111 32.0926 11.8435 8.7892 -47.6365 32.0926 113.3643 -45.6810 7.1108 36.6696 11.8435 -45.6810 324.4745 -51.8259 3.5117 8.7892 7.1108 -51.8259 50.9129

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

P2 = 
   0.7224 + 0.1105i
   0.7224 - 0.1105i
   0.6348 + 0.2625i
   0.6348 - 0.2625i
   0.3410 + 0.0000i

%Compare K2 and K1, select which produces better results for your real
%world physical system through H2 = A-B*K2
%% For just LQI the integrator and states gain:
K_x = Klqi(1:5);
K_i = Klqi(6);
dt = s_time;  % time step;
T = 100;  % sim time in [sec]
%% For just LQR: System's Hamiltonian
H = A - B*K1

H = 5×5

0.7104 0.0194 -0.2595 -0.1524 0.2670 0.1892 0.3864 -0.3014 0.4381 -0.0459 -0.0962 0.5019 0.7347 -0.0399 0.1894 0.0144 0.1372 -0.0602 0.8170 -0.0921 -0.1094 0.4205 -0.1416 -0.3972 0.3907

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

H2 = A - B*K2

H2 = 5×5

0.7118 0.0192 -0.2552 -0.1512 0.2687 0.1833 0.3874 -0.3201 0.4329 -0.0535 -0.0942 0.5016 0.7410 -0.0381 0.1919 0.0171 0.1367 -0.0517 0.8194 -0.0887 -0.1055 0.4198 -0.1293 -0.3937 0.3958

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%% Now when the optimal quadratic gain matrix K is found for your
%% embedded implementation you will need the plant(environment) under control
%% state vector X estimator (it can be Kalman or another more handy and simple,
%% which =?(additional research question m to answer if stakeholder/customer
%% hires me). Then The control law of your system under LQR control is ultimately
%% u = setpointreferenceinputsignal - K*Xest. Very impressive gain and phase margines if
%% we find as precise Xest-X -> 0 and and also for homogeneous and in-homogeneous cases.
% 
Gn = 1; %you may experiment with different gains
sys3 = ss(H,B,Gn*C,D,s_time);
isstable(sys3)

ans = logical
   1

figure
step(sys3);

stepinfo(sys3)

ans = struct with fields:
         RiseTime: 0.1400
    TransientTime: 0.4200
     SettlingTime: 0.4200
      SettlingMin: 0.1474
      SettlingMax: 0.1627
        Overshoot: 0.0013
       Undershoot: 5.9286
             Peak: 0.1627
         PeakTime: 0.8400

hasdelay(sys3)

ans = logical
   0

tdc3 = totaldelay(sys3)

tdc3 = 0

sysndc3 = absorbDelay(sys3)

sysndc3 =
 
  A = 
             x1        x2        x3        x4        x5
   x1    0.7104   0.01942   -0.2595   -0.1524     0.267
   x2    0.1892    0.3864   -0.3014    0.4381  -0.04587
   x3  -0.09623    0.5019    0.7347  -0.03993    0.1894
   x4   0.01441    0.1372  -0.06018     0.817  -0.09214
   x5   -0.1094    0.4205   -0.1416   -0.3972    0.3907
 
  B = 
             u1
   x1  -0.06913
   x2    0.2981
   x3   -0.1013
   x4   -0.1356
   x5   -0.1971
 
  C = 
             x1        x2        x3        x4        x5
   y1    0.3955   0.08747    0.1275  -0.05117  -0.04206
 
  D = 
       u1
   y1   0
 
Sample time: 0.02 seconds
Discrete-time state-space model.

sysxc3 = pade(sys3,nstates)

Warning: The model is already discrete. For discrete models, the "pade" command maps all delays to poles at z=0.

sysxc3 =
 
  A = 
             x1        x2        x3        x4        x5
   x1    0.7104   0.01942   -0.2595   -0.1524     0.267
   x2    0.1892    0.3864   -0.3014    0.4381  -0.04587
   x3  -0.09623    0.5019    0.7347  -0.03993    0.1894
   x4   0.01441    0.1372  -0.06018     0.817  -0.09214
   x5   -0.1094    0.4205   -0.1416   -0.3972    0.3907
 
  B = 
             u1
   x1  -0.06913
   x2    0.2981
   x3   -0.1013
   x4   -0.1356
   x5   -0.1971
 
  C = 
             x1        x2        x3        x4        x5
   y1    0.3955   0.08747    0.1275  -0.05117  -0.04206
 
  D = 
       u1
   y1   0
 
Sample time: 0.02 seconds
Discrete-time state-space model.

[DelayFreeSys3,tau3] = getDelayModel(sys3)

DelayFreeSys3 =
 
  A = 
             x1        x2        x3        x4        x5
   x1    0.7104   0.01942   -0.2595   -0.1524     0.267
   x2    0.1892    0.3864   -0.3014    0.4381  -0.04587
   x3  -0.09623    0.5019    0.7347  -0.03993    0.1894
   x4   0.01441    0.1372  -0.06018     0.817  -0.09214
   x5   -0.1094    0.4205   -0.1416   -0.3972    0.3907
 
  B = 
             u1
   x1  -0.06913
   x2    0.2981
   x3   -0.1013
   x4   -0.1356
   x5   -0.1971
 
  C = 
             x1        x2        x3        x4        x5
   y1    0.3955   0.08747    0.1275  -0.05117  -0.04206
 
  D = 
       u1
   y1   0
 
Sample time: 0.02 seconds
Discrete-time state-space model.

tau3 =

  0x1 empty double column vector

% Convert state-space model to transfer function
sys_tf  = tf(sys3);           
[n, d]  = tfdata(sys_tf, 'v'); 
% Pre-filter
% !Pre-filter, let's behave professionally, the original idea of such
% fitler is of Sam Chak from here and so I must cite him and credit his
% contribution: https://uk.mathworks.com/matlabcentral/profile/authors/5570092
%https://uk.mathworks.com/matlabcentral/answers/2112651-i-am-trying-to-find-the-lqr-controller-system-transfer-function-but-getting-error-as-arrays-have#answer_1449181
AdjGain = 0.02;
% with my additions of more scripting/coding control
fdst = length(d)-1

fdst = 5

fdend = length(d)

fdend = 6

fnst = 1

fnst = 1

fnend = length(n)

fnend = 6

Gf = AdjGain*tf(d(fdst:fdend),n(fnst:fnend),s_time)   

Gf =
 
                      0.01367 z - 0.001616
  -------------------------------------------------------------
  0.001044 z^4 - 0.01182 z^3 + 0.0379 z^2 - 0.05049 z + 0.02531
 
Sample time: 0.02 seconds
Discrete-time transfer function.

       
% Display the transfer function of the  Compensated System
disp('Transfer Function of the Compensated System:');

Transfer Function of the Compensated System:

Gcc = minreal(Gf*sys_tf)        

Gcc =
 
                     0.01367 z - 0.001616
  -----------------------------------------------------------
  z^5 - 3.039 z^4 + 3.713 z^3 - 2.264 z^2 + 0.6835 z - 0.0808
 
Sample time: 0.02 seconds
Discrete-time transfer function.

ms = minreal(Gcc);
zms = zero(ms)     % closed-loop zeros

zms = 0.1182

pms = pole(ms)     % closed-loop poles

pms = 
   0.6870 + 0.2420i
   0.6870 - 0.2420i
   0.6773 + 0.0000i
   0.6323 + 0.0000i
   0.3556 + 0.0000i

%Step response of the resultant system to step in amplitude 12 applied at
%time = 1 sec
% replaces loop:
tsim = 10;
setpointamp = 24.5; %Bar or Newton or mm
setpointapptime = 0.9; % sec
defaultsetpointpos = 0; % Bar or Newton or mm
Conf = RespConfig(Bias=defaultsetpointpos,Amplitude=setpointamp,Delay=setpointapptime)

Conf = 
  RespConfig with properties:

                Bias: 0
           Amplitude: 24.5000
               Delay: 0.9000
        InitialState: []
    InitialParameter: []

figure
step(Gcc,[0 tsim], Conf)
title('Plant+Controller LQR Closed Loop system step response')

stepinfo(Gcc)

ans = struct with fields:
         RiseTime: 0.1800
    TransientTime: 0.3600
     SettlingTime: 0.3600
      SettlingMin: 0.9371
      SettlingMax: 1.0072
        Overshoot: 8.2561e-04
       Undershoot: 0
             Peak: 1.0072
         PeakTime: 0.7800

hasdelay(Gcc)

ans = logical
   0

tdc = totaldelay(Gcc)

tdc = 0

sysndc = absorbDelay(Gcc)

sysndc =
 
                     0.01367 z - 0.001616
  -----------------------------------------------------------
  z^5 - 3.039 z^4 + 3.713 z^3 - 2.264 z^2 + 0.6835 z - 0.0808
 
Sample time: 0.02 seconds
Discrete-time transfer function.

sysxc = pade(Gcc,nstates)

Warning: The model is already discrete. For discrete models, the "pade" command maps all delays to poles at z=0.

sysxc =
 
                     0.01367 z - 0.001616
  -----------------------------------------------------------
  z^5 - 3.039 z^4 + 3.713 z^3 - 2.264 z^2 + 0.6835 z - 0.0808
 
Sample time: 0.02 seconds
Discrete-time transfer function.

%[DelayFreeSys1,tau1] = getDelayModel(Gf*sys_tf)

%% PID
P1   = tf(sys1)

P1 =
 
  0.001044 z^4 - 0.01182 z^3 + 0.0379 z^2 - 0.05049 z + 0.02531
  -------------------------------------------------------------
   z^5 - 3.659 z^4 + 5.401 z^3 - 4.011 z^2 + 1.491 z - 0.2201
 
Sample time: 0.02 seconds
Discrete-time transfer function.

[C, info] = pidtune(P1, 'PIDF2')

C =
 
                        Ts  
  u = Kp (b*r-y) + Ki ------ (r-y)
                        z-1 

  with Kp = 0.442, Ki = 2.92, b = 1, Ts = 0.02
 
Sample time: 0.02 seconds
Discrete-time 2-DOF PI controller in parallel form.

info = struct with fields:
                Stable: 1
    CrossoverFrequency: 2.3678
           PhaseMargin: 60.0000

C2tf = tf(C)

C2tf =
 
  From input 1 to output:
  0.4418 z - 0.3835
  -----------------
        z - 1
 
  From input 2 to output:
  -0.4418 z + 0.3835
  ------------------
        z - 1
 
Sample time: 0.02 seconds
Discrete-time transfer function.

Cr = C2tf(1)

Cr =
 
  0.4418 z - 0.3835
  -----------------
        z - 1
 
Sample time: 0.02 seconds
Discrete-time transfer function.

Cy = C2tf(2)

Cy =
 
  -0.4418 z + 0.3835
  ------------------
        z - 1
 
Sample time: 0.02 seconds
Discrete-time transfer function.

Gclp = Cr*feedback(P1,Cy,+1)

Gclp =
 
  0.0004612 z^6 - 0.006086 z^5 + 0.02691 z^4 - 0.05813 z^3 + 0.06739 z^2 - 0.04026 z + 0.009708
  ---------------------------------------------------------------------------------------------
       z^7 - 5.659 z^6 + 13.71 z^5 - 18.45 z^4 + 14.86 z^3 - 7.146 z^2 + 1.891 z - 0.2104
 
Sample time: 0.02 seconds
Discrete-time transfer function.

Gcl = minreal(Gclp);
zcl = zero(Gcl)     % closed-loop zeros

zcl = 
   7.1078 + 0.0000i
   1.3136 + 0.6463i
   1.3136 - 0.6463i
   1.5918 + 0.0000i
   1.0000 + 0.0000i
   0.8680 + 0.0000i

pcl = pole(Gcl)     % closed-loop poles

pcl = 
   0.7309 + 0.3420i
   0.7309 - 0.3420i
   1.0000 + 0.0000i
   0.9419 + 0.0965i
   0.9419 - 0.0965i
   0.9224 + 0.0000i
   0.3906 + 0.0000i

% Plot the results
subplot(3,1,1)
step(sys1,[0 tsim], Conf),       grid on, legend('Original system',    'location', 'east')
subplot(3,1,2)
step(Gcl,[0 tsim], Conf), grid on, legend('2DOFPIDF Compensated system', 'location', 'east')
subplot(3,1,3)
step(Gcc,[0 tsim], Conf), grid on, legend('LQR Controlled system', 'location', 'east')

%full, complete component so far. Clear improvement in magnitude and
%settling time response of the discrete time closed loop under the
%specified control.
% 
 pause

Warning: Pause with no arguments is not supported in non-interactive contexts.

%[K,S,e] = lqry(sys,Q,R,N)
%% More robust gaussian servo controller:
nx = 5;    %Number of states
ny = 5;    %Number of outputs
%To improve the sharpness and settling of the LQG1 step response, you can try adjusting the weighting matrices and controller design parameters. Here are some suggestions to enhance the performance of the KLQG1 controller:
    %Increase the weight on the state variables in the QXU matrix to prioritize the control effort's impact on the system states.
    %Tune the weights in the QWV matrix carefully, especially the one related to the noise covariance (Rn) to reduce its influence on the closed-loop response.
    %Experiment with different values for the QI matrix to adjust the integral control's contribution to the overall controller behavior.
    %Consider refining the controller design by incorporating additional design constraints or performance specifications tailored to your specific requirements.
%By iteratively adjusting these parameters and observing their effects on the closed-loop system response, you can fine-tune the KLQG1 controller to achieve a sharper and more settled step response. 
% Define the different values of Rn to test
Rn_values = [0.01, 0.05, 0.1, 0.5, 1];
P = sys1;
QI = 5.5 * eye(ny);
Qn = Q2; % + 50 * [0.05 0.03 0 0 0; 0.03 0.02 0 0 0; 0 0 0.01 0 0; 0 0 0 0.01 0; 0 0 0 0 0.01];
R = 0.1*sqrt(Qn);
QXU = 30*blkdiag(eye(nx), R);


%Weighting matrices Q, R, and N, which define the tradeoff between regulation performance (how fast x(t) goes to zero) and control effort.
%x = [0;0;0;0;0];
%x_i = x;
%y = [];
%u = [];
%ref = [];
%states = [];
%r=0;
%Todo: simulate it in loop:
%noise = 0.001;
%ToDo next is not yet fully full: not full, under development.
%for t = time
%    if t == 50
%        r = 5;
%    end

%    u_t = -K_x * x - K_i * x_i;

%    u_t(abs(u_t)>abs(max(C))) = max(C);
%    x(abs(x)>abs(max(r))) = r;

%    x_dot = A * x + B' * u_t;
%    xn = x + x_dot * dt;
%    y_t = C * xn + D*u_t;
    %Trapezoidal integrator (ToDo: try other discrete integrators approximations, like backward and trapezoidal)
%    x_i_dot = r - y_t + r - C * x + D*u_t;
%    x_i = x_i + x_i_dot * dt/2;
%    y = [y, y_t];
%    u = [u, u_t];
%    ref = [ref, r];
%    x = xn;
%    states = [states, x];
%end

%figure;
%subplot(3,1,1);
%plot(time, ref, 'r--', 'LineWidth', 1.5);
%hold on;
%plot(time, y, 'b', 'LineWidth', 1.5);
%xlabel('Time (miliseconds)');
%ylabel('Output y');
%title('System Output vs Reference Signal');
%legend('Reference', 'Output');
%grid on;
%subplot(3,1,2);
%plot(time, u, 'g', 'LineWidth', 1.5);
%xlabel('Time (miliseconds)');
%ylabel('Control Effort u');
%title('Control Effort');
%grid on;
%subplot(3,1,3);
%plot(time, states, 'LineWidth', 1.5);
%xlabel('Time (miliseconds)');
%ylabel('States');
%title('States over Time');
%legend('x1', 'x2', 'x3', 'x4', 'x5');
%grid on;

%Option 2, construct positive semidefinite matrices:
%Nn = 1;
%[kest,L,P2] = kalman(sys1,Qn,Rn,Nn);
%Connect the  and sys to form the working loop, it should bring more performant control:
%K4 = lqi(sys1,Q,R)
%servocontroller = lqgtrack(kest, K4);
%C=servocontroller;
%set(C,'InputName','e')
%set(C,'OutputName','u')
%set(P,'OutputName','y')
%set(P,'InputName','u')
%Sum_Error = sumblk('e = r - y');
%clsys2 = minreal(tf(connect(C,P,Sum_Error,'r','y')))

%figure
%step(clsys2)
%title('Constructed closed loop system response 2')

%Option 3:
%C = lqgtrack(kest,K4)
%C = lqgtrack(kest,K4,'2dof')
%C = lqgtrack(kest,K4,'1dof')

%Fcl = minreal(series(P, C));
%serc = minreal(feedback(Fcl,sumblk));

%Constructed from needing help prgm. by
    %https://independent.academia.edu/PMazniker
    %[phone], https://join.skype.com/invite/oXnJhbgys7oW
    %https://diag.net/u/u6r3ondjie0w0l8138bafm095b
    %https://github.com/goodengineer
    %https://orcid.org/0000-0001-8184-8166
    %https://willwork781147312.wordpress.com/portfolio/cp/
    %https://www.youtube.com/channel/UCC__7jMOAHak0MVkUFtmO-w
    %https://nanohub.org/members/130066
    %https://pangian.com/user/hiretoserve/
    %https://substack.com/profile/191772642-paul-m

(here obvious delay reduction demonstrated as well as amplification and 1 loop reduction (step with config instead of for loop))

Simplest design of LQR, LQI solution presented:

because you created discrete time state space model and the lqr with matrices deals with continuous time.

What is the value of N and S, E Ricatty matrices for more precise control(agent actions) for your specific plant(environment)?

[K,S,P] = lqr(sys,Q,R,N) for discrete systems.

[K,S,e] = lqi(SYS,Q,R,N) for discrete systems.

[K,S,P] = dlqr(A,B,Q,R,N) for discrete systems only.

If you putted there matrices (A,B,C,D) it is for continuous systems.

Also the disturbance and noise impact-cost to analyze more and the system can be more tuned (settling time, overshoot, phase and gain margins, response times may be improved).

I plan to continue to elaborate it further to most precise completion. There are a few solutions availble to this.

What do you mean by the plant? What is the real world non-linear time varying (little time varying) initial the plant model? At which point linearized? How did you obtain the A,B,C,D, matrices? Which plant is it exactly (is it servomotor+pedal+the rest of braking system or is it just braking pedal and rest of braking system or just servomotor)?

What is the main goal to control mostly the braking pedal relative displacement (and so 1-to-1 mapping the servo position as the main set point) or braking pedal force excerted to control force applied on pedal or just desired braking pressure (main set-point)?

Should I make the few versions of the real world setup plants for you?

In simulation you are on time, in reality it is discrete time (quantization and sampling impacts are added).

So on discretization.

Correction to your post:" I have neve had any issues with goal reaching".

What are the discrete controller and plant precise parameters as in your real world system?

Correction to your post:

Your wrote in post :"servomechanism LQI controller."

But in code there is LQR controller.

I added both LQR and LQI design for discrete.

Have you tried the LQG control I proposed as more fitted clear in previous your similar post?

Have you tried PID tuning?

So far you made in script. I've just started to ammend, improve your the code script, the system and controller design and preparing of Simulink diagram, block-scheme for it...

Hope to help, can you await till I complete and upload it? Here is good educational summary of available controllers for you(I have found a few more, not yet presented in the TCE Mathworks matlab self-paced curses, which is fPID, LQG,...more I can find/look up and put effort to understand if there will be found stakeholders(customers,clients,invsestors) to order the Project service and products (both digital products and processes contribution which will lead to products and services)

Here I've created novel model and sim with your A,B,C,D, matrices and step input and white noise disturbances in actuator (control agent) output(actions) and plant output (sensor noise) to model it as close to real world and with PID(z) controller with Forward Euler and Trapezoidal Integration and filtered derivative, not saturated, 3 signals(control effort, yellow - plant resultant output, setpoint) presented. As we can see the hope is robust, the plant got stabilized indeed in full with the PID and the output is even over-gained, while delay in output is still up to 1 sec even in the PID control case, which is quite large and for it others use specific sluggish systems treatment techniques, which I think to cover next. The design is modular, hypothethical models also presented in it in case to grow it modular.

[Wiktoria]

Hello!

Thank you for your answer! I have found a research that implemented LQI based on augmenting the matrices so that the integral term is taken into consideration and then just using lqr function. From what I understood, it I augment the matrices so I take the integral into consideration, it is basically LQI. Perhaps that was wrong?

I haven't tried simulink because I have no idea how to convert the code from simulink to embedded (in my case, STM32). Additionally, the servo braking is a part of a bigger project, so I am a bit confused with how to mend it all together.

My plant is basically a bolide (formula student car). I got the state space matrices after obtaining response to step from system (but step was from 0 bars to 24 bars) and then running System Identification Toolbox to estimate state space model. The plant is servomechanism connected to a brake pedal via steel wire and it gets requested brake pressure via CANbus.

The goal is to implement the controller with a good transient response (so no excessive or sharp changes in control signal) and to make it robust, reliable and fast. We have used a PID controller before but - at least from what I got to know - it didn't give satisfactory performance. My task was to implement a different one (like LQR, LQI or something else) to compare the solutions and performence between the new solution and the previous one.

I have tried to implement the LQG controller you have mentioned, but I am confused with the equation the control signal should have. I tried to compute it, looked up on MATLAB guidelines, but when I implemented it, it was unstable, so I tried to fix my current solution. I have tried to compute the u_t using gain counted by dlqr function, as + noise.

Since my goal is to implement a different controller than PID, I haven't tried PID tuning.

My state space model is for sure discrete. I have check it and MATLAB returns that it is a discrete time, sampling 0.02s.

The servomechanism works for pulse lengths from 900 us to 1300 us. In current solution, the controller computes the control signal (as position), remaps it so it's between 0 to 90 degrees and then computes the pulse length for the computed, rescaled position. The delay between the servo receiving the command and reacting to it is around 10 to 20 ms, resolution of the servo is <2us. In embedded system (STM32) servo actuation is set up on 50Hz timer (so every 0.02s). First is checks if there is a need for actuation (by checking if error between goal pressure and current pressure is acceptable) and if it isn't, the controller starts to compute the output and servo moves until it reaches the goal and is within the acceptable error range.

[Pavl M.]

OK, it is doable (both 2DOF fractionalPID and LQ-family of controllers (including LQR-PID hybrids, and loop recovery can be tried on it), when I complete I will upload it:

Loop Recovery techniques are very promissing since allow more accuracy than Kalman X state estimator. Let's check it.

Please check recent update to the initial answer here, fast and robust K gain for LQR controller and 2DOFPID and tracking found. Let it possible I will add more to it.

Similar very related question:

https://www.mathworks.com/matlabcentral/answers/530608-how-do-i-get-my-desired-settling-time-or-rise-time-using-lqr-which-things-i-need-to-change-to-get-n?s_tid=answers_rc1-2_p2_MLT

Carefull selection of Q,R,N matrices are similar in goal vs the process of PID tuning.

Let's see maybe you can use the controller (smart agent for which the good K gains finding we have completed) with your system without complex Particle Filter or Kalman state X estimator if to rely on your matrices A,B,C,D (please check are you sure about the matrices, whether they work in embedded code?) solely then the state X can be computed inside the embedded code from equation

X[n+1] = A*X[n]+B*U[n],

then the control output can be computed as:

LQR:

U = setpoint - Klqr1(1...Nstates)*X(1..Nstates) (designed novel so far (see code K1 gain) so far with no augmentation)

(most rapid & accurate, it meets your requirements and goals as for rise time, overshoot, gain and settling), but how is it noise prone & disturbance withstanding in real world harsh automotive environments?)

U = (setpoint - y_from_real_phys_system)*Klqr2(1) - (Klqr2(2).. Klqr2(N))*(X2 .. XN) or (X1..XN-1)(can be augmented in a few ways: pre-augmented, in-augmented, post-augmented and augmented from beginning and from end)

or U = F1(setpoint - y_from_real_phys_system) - (Klqr1...KlqrN)*(X1..XN), where F1 is a novel function/filter to design with better performance hope.

LQI:

u = DiscreteIntegratorApproximator(r - y_from_real_phys_system)*Klqi(N+1) - Klqi(1..N)*X(1..N)

or

u = DiscreteIntegratorApproximator(r - y_from_real_phys_system)*Klqi1 - Klqi(2..N+1)*X(1..N). (depending on augmentation from beginning or augmentation at end)

or

u = DiscrIntegratorApprox(r - (y_from_real_phys_sys+disturbance)*F1)*K1 - K2*X(preselected states acc. to states physical sense) +disturbance*F2

where F1 is some semi-stochastic function of harshnessin real environement where plant and controller act in closed loop, catched at actuator disturbance hotstop,

K1 - is part for the new state of augmented or even not augmented gain matrix computed by some most fitted LQ algorithm,

K2 - is part for selected exisiting real physical plant states of augmented or even not augmented gain matrix computed by some most fitted LQ strategy(by algebraic, matrix and calculus operations).

[K1 K2] = K cane be computed by 1 algorithm or K1 and K2 can be also calculated separatedly by different math. programming objections emending.

where F2 is some semi-stochastic function of harshness in real environement where plant& controller act in closed loop, catched at real world physcial plant sensor/transducer disturbance hotspot.

The control u computed on your embedded board relying on the actual real physical system estimation matrices A,B,C,D, hard coded in your embedded code and using actual physical system output y, the control u computed you just connect via appropriate transducer (signal conditioner) to your actual physical system actuator (servo motor or any other valid actuator).

Let's focus on this workflow.

The control u computed on your embedded board relying on the actual real physical system estimation matrices A,B,C,D, hard coded in your embedded code and using actual physical system output y, the control u computed you just connect via appropriate transducer (signal conditioner) to your actual physical system actuator (servo motor or any other valid actuator).

Let's focus on this workflow for next few monthes and we will see also can we obtain in practice better than PID control with LQR (in theory LQR is better than PID, but whether better than the proposed 2DOFPID, with best phase and gain margins), but in practice due to the X state observer additional design let's see what we get in practice in your embedded board in C,C++ controller running compiled (in mind) (please keep your the setup, I will not steal it, while already today others many stolen from me more locally, that is why I request to find material substance support investor, granter to help only financially by Money Gram, RIA or Western Union remiitance to Pavlo Mazniker) . It is very valued industry and academic research question with many unknowns to commerciallize.

OK, as for the LQ RAE goal conroller(agent), it is feasible of course, for just LQR/LQI it must needs to

design implementable online (real-time) closed loop Kalman state estimator as much accurate and robust, while the K gains we have already found, I have started Simulink diagram blockscheme design of LQ with different Kalmans to select best one and also thinking more robust and more implementable on digital hardware (Which is better MCU, CPLD, FPGA, custom ASIC? for future applications?) than Kalman.

I tried today also to close and run loop with just 1 feedback gain element of value operation_on(K) and in feedforward integrator + the operation_on(K) gain, but yet the output is going down, sign that more complex but not complicated design is on agenda.

The question is that from time to time the LQR-LQG controllers found different interpretations in academic research and whether online(real time) closed loop Kalman gain estimator is required if we already pre-computed the K gains?

We need to distinguish between Kalman state X estimator and Kalman K gain estimator

and 3 main families of Ricatti Optimization LQ Designs:

LQG with Kalman K gain vector estimator and 1 integrator

LQR wich is just I controller from PID tuned using Ricatti Optimization with balance matrices carefully selected as in the rest of the kinds. Here also experiments with vector/matrix of K gains transformation as a feedback gain or

feedforward as I controller or a combination of both can be conducted if you contact and hire me.

LQI is a controller with gain, integrator and Kalman plant internal state(X) estimation

LoopTransferRecovery (LTR) technique,

What are your heuristics to choose Q,R,N matrices?

In previous the PID closed loop Simulink simulation presented, let me correct my the PID for comparison: there should be - in sum block on output signal on feedback, that picture were produced with + sign there.

It can be doable in TCE Simulink of course, have you tried Scicos/Scilab?

Are you sure about 900 to 1300 [us = 0.000001 sec] and not [0.001 sec = ms], please check time units consistency:

please check:

900 us to 1300 us. In current solution, the controller computes the control signal (as position),

remaps it so it's between 0 to 90 degrees and then computes the pulse length for the computed,

rescaled position. The delay between the servo receiving the command and reacting to it is around 10 to 20 ms,

resolution of the servo is <2[us]

Are you sure about <2[us]? I think change your us to ms.

Is the servo reacting delay already resembled in pre-obtained the A,B,C,D matrices (are needed to implement the actual controller(agent) logic in embedded C) or should we model / regard it additionally (as *exp(... delay_value))?

Answer 3

[Paul]

Following up from this comment, we have established that the plant is modeled in disctete-time with sampling period Ts = 0.02 s.

% System identification
A = [0.748153090242338 -0.0770817542371630 -0.272770107988345 -0.0680287168856025 0.274040229219026;
    0.0264773037738115 0.802405697397900 -0.244152268701560 0.0743061296316669 -0.0764153326240192;
    -0.0409228451430832 0.360565865830120 0.715237466131514 0.0837023187613815 0.199733910116335;
    0.0884449413535669 -0.0520801836735802 -0.0862407832241391 0.982500599129848 -0.0782458817344473;
    -0.00176167556838596 0.145291619650476 -0.179472954855322 -0.156570713297809 0.410924513958431];
B = [-0.0691336950664988;
    0.298054008136356;
    -0.101282668881553;
    -0.135571075515549;
    -0.197136001368224];
C = [0.395474090810162 0.0874731850181924 0.127527899108947 -0.0511706805106379 -0.0420611557807524];
D = 0;
Ts = 0.02;
sys = ss(A, B, C, D, Ts);

However, the rest of the code in the Question, i.e., the augmentation, the call to lqr, and the simulation, is treating the model as if its in continuous-time, so that all needs to be corrected.

It appears the objective is to use the LQR technique with integral control, for which we have two options. Option 1 is to manually augment the plant with a discrete-time integrator and then call lqr with the augmented LTI object (not Aaug and Baug), Option 2 is to use lqi and let the software do the work for us.

Let's use option 2.

Q = 70*eye(size(A,1)+1);
R = 5;
K = lqi(sys,Q,R);

Now that we have our gains, we have to build the closed loop system based on those gains. We have some options here. Option 1 is to manually manipulate the models to come up with a closed-form expression. Option 2 is to use connect so we can more intuitively represent the system as a block diagram, and it also provides more flexibility going forward when we want to do other things with the model, like add models for a sensor or specficy breakpoints for stability margins. Let's use option 2

The plant model presently has only one output, y, but we also need to output the state variable for state feedback

sys = ss(A,B,[C;eye(5)],0,Ts,'InputName','u','OutputName',["y", "x("+(1:5)+")"]);

Define the sum block that forms the error signal

errsum = sumblk('e = r - y');

Define the discrete time integerator using the equation shown at lqi

intsys = ss(1,Ts,1,0,Ts,'InputName','e','OutputName','xi');

Define the gains

K_x = ss(K(1:5),'InputName','x','OutputName','ux');
K_i = ss(K(6),'InputName','xi','OutputName','ui');

Sum of the control signals from the state feedback and the integral control

usum = sumblk('u = -ui - ux');

I think you're going to want a stand alone model of the compensator, so let's build that first

syscomp = connect(errsum,intsys,K_x,K_i,usum,["r";"y";K_x.InputName],'u');

Put it all together

syscl = connect(sys,syscomp,'r','y');

Now, for simulation we can either write our own loop, or we can just use step. Compare the reponse of the plant to the closed-loop system

figure
step(sys('y','u'),syscl)

The slow response of the closed loop system is due to the selection of Q and R. I don't know how you came up with those or what your design goals are, so won't suggest how to modify those. I can give you an idea of one approach that might be worth pursuing if you'd like.

[Wiktoria]

Thank you so much for your help! I can fully see now why my approach wasn't good. My main concern though is that I need to develop the loop by myself, since the code will be rewritten to C. I wanted to create a loop step by step, like it is done in my (now wrongly implemented) code. Additionally, the Q and R matrices are choosen incorrectly - they used to work better with my previous solution, but now they require changing. I plan on using the Bryson's rule to determine those. My goal is to implement a fast and reliable system without much change in the control signal.

I will be thankful for any ideas and I'll gladly learn about a new approach I can use. Thank you!

[Paul]

For selecting Q and R, Bryson's rule might be a good starting point if you have a good understanding of the physical interpreration of the state variables, I suppose.

Here's another approach you can experiment with.

% System identification
A = [0.748153090242338 -0.0770817542371630 -0.272770107988345 -0.0680287168856025 0.274040229219026;
    0.0264773037738115 0.802405697397900 -0.244152268701560 0.0743061296316669 -0.0764153326240192;
    -0.0409228451430832 0.360565865830120 0.715237466131514 0.0837023187613815 0.199733910116335;
    0.0884449413535669 -0.0520801836735802 -0.0862407832241391 0.982500599129848 -0.0782458817344473;
    -0.00176167556838596 0.145291619650476 -0.179472954855322 -0.156570713297809 0.410924513958431];
B = [-0.0691336950664988;
    0.298054008136356;
    -0.101282668881553;
    -0.135571075515549;
    -0.197136001368224];
C = [0.395474090810162 0.0874731850181924 0.127527899108947 -0.0511706805106379 -0.0420611557807524];
D = 0;
Ts = 0.02;
sys = ss(A, B, C, D, Ts);
zpk(sys)

ans =
 
     0.0010438 (z-7.108) (z-1.592) (z^2 - 2.627z + 2.143)
  ----------------------------------------------------------
  (z-0.4701) (z^2 - 1.708z + 0.7385) (z^2 - 1.482z + 0.6339)
 
Sample time: 0.02 seconds
Discrete-time zero/pole/gain model.

Define the quadratic cost to be minized as: H = a*y^2 + b*yi^2, where y is the output of the plant and yi is integral of y, whics is the integrator state in the control, xi. Then, using the fact that y = C*x, we have

H = a*x' * C' * C * x + b*xi^2,

which can be combined as

a = 10; b = 500;
Q = blkdiag(a*C'*C,b);
R = 1;  % everything scales with R, so just set it to one
K = lqi(sys,Q,R);

sys = ss(A,B,[C;eye(5)],0,Ts,'InputName','u','OutputName',["y", "x("+(1:5)+")"]);
errsum = sumblk('e = r - y');
intsys = ss(1,Ts,1,0,Ts,'InputName','e','OutputName','xi');
K_x = ss(K(1:5),'InputName','x','OutputName','ux');
K_i = ss(K(6),'InputName','xi','OutputName','ui');
usum = sumblk('u = -ui - ux');
syscomp = connect(errsum,intsys,K_x,K_i,usum,["r";"y";K_x.InputName],'u');
syscl = connect(sys,syscomp,'r','y');
figure
step(sys('y','u'),syscl)

I just played around a bit with a and b to get a reasonable looking response. Not sure what you're looking for, but with this approach there are only two parameters to iterate over to try to meet the design objectives.

There are probably other methods for selecting Q and R.

As for simulating everything in a loop, you have the difference equations of the controller and the plant as defined by their state space models that you can iterate. But, from what I understand of your project you really only need to implement syscomp in C code.

Also, I'm wondering about the plant model. It kind of looks like the plant has about 10 cycles of delay, which is somehow captured in the A/B/C/D dynamics. Maybe the system identification should try capture the delays explicitly. As it stands now, it's possible that you won't have all of the states of the identified system available for feedback, but I don't really know what your system is or what your modeling, so my concerns could be unfounded.

[Wiktoria]

To be honest, I actually need the loop written step by step like in the code I provided, since it makes it easier for me to understand as well as I have no idea how to switch from MATLAB to embedded C when I have so many functions used. Like I mentioned, I am a beginner and although I find control theory fascinating, it sure is a lot to unfold.

In reality just a subsystem supposed to control the servo - all I receive is the goal pressure I need to achieve and I have to control the servomechanism accordingly. Goal pressure is set via CAN. I have heard that there might be a small delay caused by both the information transfer and servo reaction, but heard not to worry as of now. I will think about trying to capture the delays in model explicitly though, I believe I can do that in System Identification Toolbox.

Thank you for all your help! I am really grateful for the explanation you provided.

[Mathieu NOE]

@Wiktoria

Hi Wiktoria,

could you share the data for the system identification

we need both the stimulus (input ) and output

I still wonder how we get the 200 ms delay

have a great day

[Pavl M.]

No free lunch,

From where they took next quadratic cost:

H = a*y^2 + b*yi^2

( I found so far b ↑and a ↓ and r ↑and a~= -b/3, a~=b values on my paper for minimization of H).

Let me know which to do for which reward?

How it correlates with J and Ricatti initial LQ control definitions.Mainly to select the cost matrices and understand where the delay in your plant(process under control,environment of the braking path) comes, whether conversions or plant subsystem introduced it or the problem in controller?

By the agent(controller) to reduce delay by reducing t_rise(rise time) and settling_time is already not bad, and not to loose in gain, while to reduce control effort how phase and gain margins performance characteristics resemble the latency reduction/quantity?

They used before zero-poles map, root locus and bode plot to adjust the weighting matrices, mainly state cost Q each cell separately as per physical sense on the X(1:i) states, cross matrix N, control cost R state-control, interstate, noise covariance, disturbance rejection vs setpoint tracking

To compute the LQG regulator, lqg uses the commands lqr and kalman. To compute the servo-controller, lqg uses the commands lqi and kalman.

LQG is good, because it produces full controller (both LQ and state estimation) in 1 just 1 completed block and takes into considerations also the plant state noise. But so far i just found it running suboptimally. For more optimal performance there will be needing more experiments with QXU QWV QI matrices.

Another good implementation roadmap is that for simplicity of realization on embedded board you may not be needing Kalman state observer without loss of control efficicency? In so case the PID and LQR with plant states estimated from the X of your system identification efforts A B C D matrices are the best fit.

You need to take into account also usually very harsh conditions of automotive systems in real life

and so MPC with robust tracker will be needed to SOA.

Right from the start in order to alleviate painfull trial-and-error iterative design procedure, a few pre-developed, pre-trained heuristics can be used:

a. Understand your control effort energy demands/limitations, your physical system model natural sens, especially the states X sense;

b. accurate1-to-1 connecting expressions (mapping morphism) between output characteristics (such as rise time, settling time, overshoot, undershoot, stability margins, gain margins, frequency response ,damping, disturbance rejection power, control effort power) and controller design parameters;

c. Closed loop architectural constraints;

If you need more specific R&D on this domain contact me more in private.

To find also: det(sI − (A − BK)) would be good to find how rise time can be expressed by Q,R,N, K is

a novel research path(topic, theme).

Please see update to the answer, whether next LQG1 and LQG2 controllers (The proposed controllers are still yet simple!) with the QXU, QWV, QI matrices selection process is similar to Q,R,N matrices (already includes in it Kalman state estimator) produce less delay faster and more optimal settlement?:

%% More robust gaussian servo controller:
% Define the different values of Rn to test
Rn_values = [0.01, 1];
P = sys1;
% Loop through each value of Rn
for i = 1:length(Rn_values)
    
    Rn = Rn_values(i);
    
    % Tune the weighting matrices
    Qn = Q2 + 50 * [0.05 0.03 0 0 0; 0.03 0.02 0 0 0; 0 0 0.01 0 0; 0 0 0 0.01 0; 0 0 0 0 0.01];
    R = sqrt(Qn);
    QXU = blkdiag(eye(nx), R);
    QWV = blkdiag(Qn, Rn);
    QI = 0.5;
    
    % Redesign the KLQG1 and KLQG2 controllers
    [KLQG1,info1] = lqg(sys1, QXU(1:A_size_row+1, 1:A_size_column+1), QWV, QI, '1dof');
    [KLQG2,info2] = lqg(sys1, QXU(1:A_size_row+1, 1:A_size_column+1), QWV, QI, '2dof');
    
    % Construct the closed-loop systems
    C = KLQG1;
    set(C, 'InputName', 'e');
    set(C, 'OutputName', 'u');
    set(P, 'OutputName', 'y');
    set(P, 'InputName', 'u');
    Sum_Error = sumblk('e = r - y');
    clsys = tf(connect(C, P, Sum_Error, 'r', 'y'));
    clsys_m = minreal(clsys);
    
    C = KLQG2;
    set(C, 'InputName', 'e');
    set(C, 'OutputName', 'u');
    set(P, 'OutputName', 'y');
    set(P, 'InputName', 'u');
    Sum_Error2 = sumblk('e = r - y');
    clsys2 = tf(connect(C, P, Sum_Error2, 'r', 'y'));
    clsys_m2 = minreal(clsys2);
    
    % Plot the results
    figure;
    subplot(2, 1, 1);
    step(clsys_m, [0 tsim], Conf);
    title(['KLQG1 Closed-Loop System Response (Rn = ', num2str(Rn), ')']);
    
    subplot(2, 1, 2);
    step(clsys_m2, [0 tsim], Conf);
    title(['KLQG2 Closed-Loop System Response (Rn = ', num2str(Rn), ')']);
    
    % Pause to allow the plots to be viewed
   
end

(simpler controllers)

Will next function help you to estimate more accurately the state space model of your plant/environment/process to control:

https://www.mathworks.com/help/ident/ref/ssest.html

ssest estimates discrete and continues state space model

Then you can get more data for your A,B,C,D matrices, ushc us nx = 5 or 4 or more.

What compare, predict and sim TCE functions produced for you identified system?

Which is optimal nx (number of states) for your the environment/plant/process model?

Use this syntax especially when you want to estimate a state-space model using frequency-domain or frequency-response data, or when you want to take advantage of the additional information, such as intersample behavior, data sample time, or experiment labeling, that data objects provide.

Non-linearities can be modeled as:

It is quite very valuable and interesting question. This control is centralized, but when we think as Ph.D. taking into account recent vehicular systems Zonal Controllers (ZCU) how to control it indeed. Also we can think more about the decentralized control: see https://en.wikipedia.org/wiki/Witsenhausen%27s_counterexample.

To be most precise the complex system of servomechanism(which is itself actuator with encoder, driver and controller)(mostly linear) then brake pedal path (that wire and brake pedal flange)-> braking cyllinder,pneumo-hydro br. tubes, lines, hoses -> brake disks/drums/pads on wheel hub(most of non-linearities) has non-linearities in it (while if it is correctly reliably constructed it is time invariant) which particular states we can clearly distinguish in it so that it is needed to think also whether QQR, PQR will be required or can we use MPC for necessary wide range of wheel forces and braking system reaction on evolving A,B,C,D linearized system matrices (each time change linearization region and reapeating the LQ(I,R or G) control)? What is the region where the actual physical system is linearized?

[Pavl M.]

clc
clear all
close all

s_time = 0.02; %[sec], SI by default in TCE Matlab

Np = 0.01; %noise power
nstates = 5;
ninputs = 1;
ncontrols = 1;
noutputs = 1;

% Initial(original) datum posted to create agent-environment closed loop:

A = [0.748153090242338 -0.0770817542371630 -0.272770107988345 -0.0680287168856025 0.274040229219026;  
0.0264773037738115 0.802405697397900 -0.244152268701560 0.0743061296316669 -0.0764153326240192;
-0.0409228451430832 0.360565865830120 0.715237466131514 0.0837023187613815 0.199733910116335;
0.0884449413535669 -0.0520801836735802 -0.0862407832241391 0.982500599129848 -0.0782458817344473;
-0.00176167556838596 0.145291619650476 -0.179472954855322 -0.156570713297809 0.410924513958431];

B = [-0.0691336950664988;   
0.298054008136356;
-0.101282668881553;
-0.135571075515549;
-0.197136001368224];

C = [0.395474090810162 0.0874731850181924 0.127527899108947 -0.0511706805106379 -0.0420611557807524];
D = 0;

sys1 = ss(A, B, C, D, s_time);
[A_size_row, A_size_column] = size(A);

hasdelay(sys1)

ans = logical
   0

td = totaldelay(sys1)

td = 0

sysnd = absorbDelay(sys1)

sysnd =
 
  A = 
              x1         x2         x3         x4         x5
   x1     0.7482   -0.07708    -0.2728   -0.06803      0.274
   x2    0.02648     0.8024    -0.2442    0.07431   -0.07642
   x3   -0.04092     0.3606     0.7152     0.0837     0.1997
   x4    0.08844   -0.05208   -0.08624     0.9825   -0.07825
   x5  -0.001762     0.1453    -0.1795    -0.1566     0.4109
 
  B = 
             u1
   x1  -0.06913
   x2    0.2981
   x3   -0.1013
   x4   -0.1356
   x5   -0.1971
 
  C = 
             x1        x2        x3        x4        x5
   y1    0.3955   0.08747    0.1275  -0.05117  -0.04206
 
  D = 
       u1
   y1   0
 
Sample time: 0.02 seconds
Discrete-time state-space model.

sysx = pade(sys1,nstates)

Warning: The model is already discrete. For discrete models, the "pade" command maps all delays to poles at z=0.

sysx =
 
  A = 
              x1         x2         x3         x4         x5
   x1     0.7482   -0.07708    -0.2728   -0.06803      0.274
   x2    0.02648     0.8024    -0.2442    0.07431   -0.07642
   x3   -0.04092     0.3606     0.7152     0.0837     0.1997
   x4    0.08844   -0.05208   -0.08624     0.9825   -0.07825
   x5  -0.001762     0.1453    -0.1795    -0.1566     0.4109
 
  B = 
             u1
   x1  -0.06913
   x2    0.2981
   x3   -0.1013
   x4   -0.1356
   x5   -0.1971
 
  C = 
             x1        x2        x3        x4        x5
   y1    0.3955   0.08747    0.1275  -0.05117  -0.04206
 
  D = 
       u1
   y1   0
 
Sample time: 0.02 seconds
Discrete-time state-space model.

[DelayFreeSys,tau] = getDelayModel(sys1)

DelayFreeSys =
 
  A = 
              x1         x2         x3         x4         x5
   x1     0.7482   -0.07708    -0.2728   -0.06803      0.274
   x2    0.02648     0.8024    -0.2442    0.07431   -0.07642
   x3   -0.04092     0.3606     0.7152     0.0837     0.1997
   x4    0.08844   -0.05208   -0.08624     0.9825   -0.07825
   x5  -0.001762     0.1453    -0.1795    -0.1566     0.4109
 
  B = 
             u1
   x1  -0.06913
   x2    0.2981
   x3   -0.1013
   x4   -0.1356
   x5   -0.1971
 
  C = 
             x1        x2        x3        x4        x5
   y1    0.3955   0.08747    0.1275  -0.05117  -0.04206
 
  D = 
       u1
   y1   0
 
Sample time: 0.02 seconds
Discrete-time state-space model.

tau =

  0x1 empty double column vector

tau

tau =

  0x1 empty double column vector

figure
w = logspace(2,6,1000);
bode(sys1,w),grid

figure
rlocus(sys1)

figure
nyquist(sys1)

figure
nichols(sys1)

Q = 7*eye(A_size_row+1)

Q = 6×6

7 0 0 0 0 0 0 7 0 0 0 0 0 0 7 0 0 0 0 0 0 7 0 0 0 0 0 0 7 0 0 0 0 0 0 7

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Q2 = 50*eye(nstates) % Optimal Q matrix

Q2 = 5×5

50 0 0 0 0 0 50 0 0 0 0 0 50 0 0 0 0 0 50 0 0 0 0 0 50

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

R = 0.002*sqrt(Q2) % R matrix

R = 5×5

0.0141 0 0 0 0 0 0.0141 0 0 0 0 0 0.0141 0 0 0 0 0 0.0141 0 0 0 0 0 0.0141

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

QXU = Q %30*blkdiag(eye(nstates), R); % Weighting matrix QXU

QXU = 6×6

7 0 0 0 0 0 0 7 0 0 0 0 0 0 7 0 0 0 0 0 0 7 0 0 0 0 0 0 7 0 0 0 0 0 0 7

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% Loop through each value of Rn
Rn = 0.0005 % Optimal Rn value

Rn = 5.0000e-04

Qn = 344*Q2 

Qn = 5×5

17200 0 0 0 0 0 17200 0 0 0 0 0 17200 0 0 0 0 0 17200 0 0 0 0 0 17200

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

QWV = 390*blkdiag(Qn, Rn)/1e+06 % Weighting matrix QWV

QWV = 6×6

6.7080 0 0 0 0 0 0 6.7080 0 0 0 0 0 0 6.7080 0 0 0 0 0 0 6.7080 0 0 0 0 0 0 6.7080 0 0 0 0 0 0 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

QI = 310*eye(noutputs) % Weighting matrix QI

QI = 310

% Redesign the KLQG1 controller with optimal values
KLQG1 = lqg(sys1, QXU, QWV, QI, '1dof');
% Construct the closed-loop system
C = KLQG1;
set(C, 'InputName', 'e');
set(C, 'OutputName', 'u');
set(sys1, 'OutputName', 'y');
set(sys1, 'InputName', 'u');
Sum_Error = sumblk('e = r - y');
clsys = tf(connect(C, sys1, Sum_Error, 'r', 'y'));
clsys_m = minreal(clsys);
% Plot the results
tsim = 10;
setpointamp = 24; % Setpoint amplitude
setpointapptime = 1; % Setpoint application time
defaultsetpointpos = 0; % Default setpoint position
Conf = RespConfig(Bias=defaultsetpointpos,Amplitude=setpointamp,Delay=setpointapptime);
figure;
step(clsys_m, [0 tsim], Conf);
title('Optimized LQG1 Closed-Loop System Response')

stepinfo(clsys_m)

ans = struct with fields:
         RiseTime: 1.6800
    TransientTime: 3.2000
     SettlingTime: 3.2000
      SettlingMin: 0.9006
      SettlingMax: 0.9986
        Overshoot: 0
       Undershoot: 0.9477
             Peak: 0.9986
         PeakTime: 5.1800

hasdelay(clsys_m)

ans = logical
   0

td = totaldelay(clsys_m)

td = 0

sysnd = absorbDelay(clsys_m)

sysnd =
 
  From input "r" to output "y":
     0.0009387 z^9 - 0.01411 z^8 + 0.07867 z^7 - 0.2344 z^6 + 0.4249 z^5 - 0.4931 z^4 + 0.3691 z^3 - 0.1725 z^2 + 0.04578 z - 0.005247
  ---------------------------------------------------------------------------------------------------------------------------------------
  z^11 - 6.95 z^10 + 21.8 z^9 - 40.64 z^8 + 49.87 z^7 - 42.07 z^6 + 24.71 z^5 - 9.967 z^4 + 2.642 z^3 - 0.4153 z^2 + 0.02936 z - 2.86e-09
 
Sample time: 0.02 seconds
Discrete-time transfer function.

sysx = pade(clsys_m,nstates)

Warning: The model is already discrete. For discrete models, the "pade" command maps all delays to poles at z=0.

sysx =
 
  From input "r" to output "y":
     0.0009387 z^9 - 0.01411 z^8 + 0.07867 z^7 - 0.2344 z^6 + 0.4249 z^5 - 0.4931 z^4 + 0.3691 z^3 - 0.1725 z^2 + 0.04578 z - 0.005247
  ---------------------------------------------------------------------------------------------------------------------------------------
  z^11 - 6.95 z^10 + 21.8 z^9 - 40.64 z^8 + 49.87 z^7 - 42.07 z^6 + 24.71 z^5 - 9.967 z^4 + 2.642 z^3 - 0.4153 z^2 + 0.02936 z - 2.86e-09
 
Sample time: 0.02 seconds
Discrete-time transfer function.

%[DelayFreeSys,tau] = getDelayModel(clsys_m)

figure
w = logspace(2,6,1000);
bode(clsys_m,w),grid

figure
rlocus(clsys_m)

figure
nyquist(clsys_m)

figure
nichols(clsys_m)

[Wiktoria]

@Mathieu NOE I just added data based on which I did the identification. Input starts from the first moment the jump started to the stabilization. Impulse was pressure from 0 to 24,5. Thank you for your time!

@Pavl M. Thank you for your input, I will analyze all the solutions you provided with great care!

[Pavl M.]

For faster & more accurate controller(agent) design 1st it is needed more substantial investment to me to utilize less latent hardware + from where I work.

So far I found the state feedback LQG controller gain for 5 and 7 states:

And also system matrices as per yr plant:

SYS2.a =

x1 x2 x3 x4 x5 x6 x7

x1 0.895 0.2623 -0.1866 0.1984 0.07844 0.108 0.05552

x2 -0.03094 0.8867 0.4076 -0.1036 -0.2934 -0.01719 -0.3774

x3 -0.0003021 -0.04986 0.8409 0.5682 -0.02406 0.2522 -0.1313

x4 0.0008012 -0.001186 -0.1231 0.7743 -0.6385 -0.02246 -0.7087

x5 -0.0007357 -0.001279 -0.02703 0.1887 0.7693 -0.6507 0.2365

x6 0.0005011 0.002218 0.009968 -0.02466 0.4015 0.4612 -0.8678

x7 -0.0003302 -0.004122 -0.01945 0.1069 0.08493 0.584 0.4143

SYS2.b =

u1

x1 0.09945

x2 0.1125

x3 0.03447

x4 0.0157

x5 -0.007276

x6 0.001673

x7 -0.0003104

SYS2.c =

x1 x2 x3 x4 x5 x6 x7

y1 0.4477 -0.4383 0.4409 -0.3602 -0.2415 -0.1787 -0.2404

SYS2.d =

u1

y1 0.00239

Sampling time: 0.01 s

Discrete-time model.

X02 =

-7.4870

-12.6791

-13.1222

-2.7618

6.1118

-1.1226

4.5275

INFO2 =

scalar structure containing the fields:

K =

4.013654

1.682882

0.092401

0.039642

0.070633

-0.109152

0.085817

Q =

9.5378e-02 5.3092e-02 1.2289e-03 -4.9161e-03 3.9237e-03 9.1983e-05 -2.7546e-04

5.3092e-02 2.9619e-02 6.9732e-04 -2.6344e-03 2.1297e-03 1.3368e-04 -1.5972e-04

1.2289e-03 6.9732e-04 3.4250e-04 1.1667e-04 -3.9732e-04 2.3722e-04 -4.6779e-04

-4.9161e-03 -2.6344e-03 1.1667e-04 4.9481e-04 -5.0088e-04 2.2905e-04 -2.2440e-04

3.9237e-03 2.1297e-03 -3.9732e-04 -5.0088e-04 8.0060e-04 -3.6649e-04 6.1913e-04

9.1983e-05 1.3368e-04 2.3722e-04 2.2905e-04 -3.6649e-04 2.7128e-04 -3.1223e-04

-2.7546e-04 -1.5972e-04 -4.6779e-04 -2.2440e-04 6.1913e-04 -3.1223e-04 6.7760e-04

Ry = 4.7048e-03

S =

2.0891e-02

1.1679e-02

1.0796e-04

-1.0833e-03

1.0499e-03

-3.4238e-05

1.7980e-04

L = 0.068592

warning: ident: IB01AD: 5: the least squares problem with coefficient matrix r_1 [6], used for computing the weighted oblique projection (for METH = 'N'), has a rank-deficient coefficient matrix

warning: called from

__slicot_identification__ at line 178 column 35

n4sid at line 240 column 21

FormulaBolidC at line 127 column 20

SYS2.a =

x1 x2 x3 x4 x5

x1 0.9031 -0.3338 -0.09819 0.3561 -0.1299

x2 0.02733 0.9094 -0.6554 -0.255 -0.5373

x3 -0.002235 0.01998 0.9463 1.044 -0.12

x4 -0.0002459 0.003027 -0.08817 0.923 0.783

x5 0.0002854 -0.001687 0.01295 -0.1593 -0.05511

SYS2.b =

u1

x1 0.07068

x2 -0.0834

x3 0.02667

x4 -0.0001719

x5 0.000417

SYS2.c =

x1 x2 x3 x4 x5

y1 0.4478 0.494 0.4581 -0.3563 0.4465

SYS2.d =

u1

y1 0.001029

Sampling time: 0.01 s

Discrete-time model.

X02 =

-2.5318

10.3018

-0.3841

2.3940

-7.4291

INFO2 =

scalar structure containing the fields:

K =

6.427595

-2.790762

-0.561979

-0.105093

0.038435

Q =

2.2884e-01 -9.4423e-02 -1.3518e-02 -1.0758e-03 1.2826e-03

-9.4423e-02 3.9114e-02 5.3428e-03 1.6368e-04 -8.7930e-04

-1.3518e-02 5.3428e-03 1.1925e-03 5.0375e-04 5.2112e-04

-1.0758e-03 1.6368e-04 5.0375e-04 5.3102e-04 6.2056e-04

1.2826e-03 -8.7930e-04 5.2112e-04 6.2056e-04 1.0509e-03

Ry = 6.9881e-03

S =

3.9964e-02

-1.6482e-02

-2.3897e-03

-2.0688e-04

1.7263e-04

L = 0.083595

There is also rank defficiency and negtive eigenvalues somewhere in your plant model!

K1 = 1×5

-0.5461 1.3958 0.1922 -1.2207 -0.1025

with just 1 filter in place closer to where LQI integrator seats:

clc
clear all
close all

%™®©

s_time = 0.02; %[sec], SI by default in TCE Matlab
Np = 0.01; %noise power
nstates = 5;
ninputs = 1;
ncontrols = 1;
noutputs = 1;

% System identification (Replace with more accurate results having for
% example sys = ssest(u,y,nx,'Ts',Ts)

A = [0.9308795641254706	-0.3081576947452329	0.07681919308097056	-0.2735999415074268	0.1415048387781206
0.04794879453274773	0.9004446378055914	0.6453239529919812	-0.002721770386339298	-0.7354398587464541
0.00680201856642904	-0.07221587634431316	0.849665637639694	0.7817090437712618	-0.4239017527219821
0.001493977508327948	-0.005394335758601366	-0.06559371453965469	0.6336472190608932	-1.242403610676959
0.0003039848966023925	-0.000717058626562655	-0.01521374071102382	0.282880500708408	0.7552836071035043];
B = [-0.05787642519196459;
0.09313927371668497;
0.06639892459449508;
0.03679743859098843;
-0.002879342503987405];
C = [-0.3897058481184756	-0.4851455782421755	0.5154800431360154	-0.4197660633133968	-0.2795058515639982];
D = 0.0001832197139953438;

sys1 = ss(A, B, C, D, s_time);

disp('Whether the system is stable, minimum phase and proper:')

Whether the system is stable, minimum phase and proper:

isstable(sys1)

ans = logical
   1

isminphase(tfdata(tf(sys1),'v'))

ans = logical
   0

isproper(sys1)

ans = logical
   1

isallpass(tfdata(tf(sys1),'v'))

ans = logical
   0

stepinfo(sys1)

ans = struct with fields:
         RiseTime: 0.1800
    TransientTime: 1.1600
     SettlingTime: 1.1600
      SettlingMin: 0.9160
      SettlingMax: 1.0685
        Overshoot: 2.6523
       Undershoot: 0.2750
             Peak: 1.0685
         PeakTime: 0.4200

figure
w = logspace(2,6,1000);
bode(sys1,w),grid

figure
rlocus(sys1)

figure
nyquist(sys1)

figure
nichols(sys1)

hasdelay(sys1)

ans = logical
   0

td = totaldelay(sys1)

td = 0

sysnd = absorbDelay(sys1)

sysnd =
 
  A = 
               x1          x2          x3          x4          x5
   x1      0.9309     -0.3082     0.07682     -0.2736      0.1415
   x2     0.04795      0.9004      0.6453   -0.002722     -0.7354
   x3    0.006802    -0.07222      0.8497      0.7817     -0.4239
   x4    0.001494   -0.005394    -0.06559      0.6336      -1.242
   x5    0.000304  -0.0007171    -0.01521      0.2829      0.7553
 
  B = 
              u1
   x1   -0.05788
   x2    0.09314
   x3     0.0664
   x4     0.0368
   x5  -0.002879
 
  C = 
            x1       x2       x3       x4       x5
   y1  -0.3897  -0.4851   0.5155  -0.4198  -0.2795
 
  D = 
              u1
   y1  0.0001832
 
Sample time: 0.02 seconds
Discrete-time state-space model.

sysx = pade(sys1,nstates)

Warning: The model is already discrete. For discrete models, the "pade" command maps all delays to poles at z=0.

sysx =
 
  A = 
               x1          x2          x3          x4          x5
   x1      0.9309     -0.3082     0.07682     -0.2736      0.1415
   x2     0.04795      0.9004      0.6453   -0.002722     -0.7354
   x3    0.006802    -0.07222      0.8497      0.7817     -0.4239
   x4    0.001494   -0.005394    -0.06559      0.6336      -1.242
   x5    0.000304  -0.0007171    -0.01521      0.2829      0.7553
 
  B = 
              u1
   x1   -0.05788
   x2    0.09314
   x3     0.0664
   x4     0.0368
   x5  -0.002879
 
  C = 
            x1       x2       x3       x4       x5
   y1  -0.3897  -0.4851   0.5155  -0.4198  -0.2795
 
  D = 
              u1
   y1  0.0001832
 
Sample time: 0.02 seconds
Discrete-time state-space model.

[DelayFreeSys,tau] = getDelayModel(sys1)

DelayFreeSys =
 
  A = 
               x1          x2          x3          x4          x5
   x1      0.9309     -0.3082     0.07682     -0.2736      0.1415
   x2     0.04795      0.9004      0.6453   -0.002722     -0.7354
   x3    0.006802    -0.07222      0.8497      0.7817     -0.4239
   x4    0.001494   -0.005394    -0.06559      0.6336      -1.242
   x5    0.000304  -0.0007171    -0.01521      0.2829      0.7553
 
  B = 
              u1
   x1   -0.05788
   x2    0.09314
   x3     0.0664
   x4     0.0368
   x5  -0.002879
 
  C = 
            x1       x2       x3       x4       x5
   y1  -0.3897  -0.4851   0.5155  -0.4198  -0.2795
 
  D = 
              u1
   y1  0.0001832
 
Sample time: 0.02 seconds
Discrete-time state-space model.

tau =

  0x1 empty double column vector

tau

tau =

  0x1 empty double column vector

Qc= ctrb(A,B)

Qc = 5×5

-0.0579 -0.0880 -0.1196 -0.1519 -0.1851 0.0931 0.1260 0.1551 0.1727 0.1720 0.0664 0.0793 0.0718 0.0439 0.0062 0.0368 0.0219 -0.0010 -0.0191 -0.0239 -0.0029 0.0071 0.0103 0.0062 -0.0015

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

controllab = rank(Qc)

controllab = 5

x = 1;
disp('The pair (A,B) must be stabilizable')

The pair (A,B) must be stabilizable

ic = [0 0 0 0 0];

figure
step(sys1);

[A_size_row, A_size_column] = size(A);
Q = 70*eye(A_size_row+1)

Q = 6×6

70 0 0 0 0 0 0 70 0 0 0 0 0 0 70 0 0 0 0 0 0 70 0 0 0 0 0 0 70 0 0 0 0 0 0 70

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

R = 7 %*eye(ninputs);

R = 7

R2 = 5 %*eye(ninputs+1);

R2 = 5

N = ones(nstates,1)

N = 5×1

1 1 1 1 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

N2 = ones(nstates+1,1)

N2 = 6×1

1 1 1 1 1 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%S
%E
%P %matrices to find furthermore.
Q2 = 30*eye(A_size_row) % To be selected using Adaptive Dynamic Programming or 

Q2 = 5×5

30 0 0 0 0 0 30 0 0 0 0 0 30 0 0 0 0 0 30 0 0 0 0 0 30

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%similar heuristic for rubustness
%R = 5;
%S = eye(nstates,noutputs);
%E = eye(nstates,nstates);
%[~,p] = chol(R)
%% For Just LQR:
[K1,S1,P1] = lqr(sys1,Q2,R)

K1 = 1×5

-0.8107 1.8887 2.7256 2.0935 -5.7534

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

S1 = 5×5

1.0e+03 * 0.1253 -0.0671 0.0226 -0.0018 -0.0416 -0.0671 0.1621 0.0561 0.0164 0.0072 0.0226 0.0561 0.2334 0.0795 -0.3008 -0.0018 0.0164 0.0795 0.2854 -0.0978 -0.0416 0.0072 -0.3008 -0.0978 1.6786

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

P1 = 
   0.6509 + 0.6264i
   0.6509 - 0.6264i
   0.7676 + 0.2623i
   0.7676 - 0.2623i
   0.7355 + 0.0000i

T = 100;  % sim time in [sec]
%% For just LQR: System's Hamiltonian
H = A - B*K1

H = 5×5

0.8840 -0.1988 0.2346 -0.1524 -0.1915 0.1235 0.7245 0.3915 -0.1977 -0.1996 0.0606 -0.1976 0.6687 0.6427 -0.0419 0.0313 -0.0749 -0.1659 0.5566 -1.0307 -0.0020 0.0047 -0.0074 0.2889 0.7387

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% 
sys3 = ss(H,B,C,D,s_time);
isstable(sys3)

ans = logical
   1

stepinfo(sys3)

ans = struct with fields:
         RiseTime: 0.0800
    TransientTime: 0.5200
     SettlingTime: 0.5200
      SettlingMin: 0.1867
      SettlingMax: 0.2159
        Overshoot: 5.1698
       Undershoot: 1.3945
             Peak: 0.2159
         PeakTime: 0.2600

figure
w = logspace(2,6,1000);
bode(sys3,w),grid

figure
rlocus(sys3)

figure
nyquist(sys3)

figure
nichols(sys3)

hasdelay(sys3)

ans = logical
   0

td = totaldelay(sys3)

td = 0

sysnd = absorbDelay(sys3)

sysnd =
 
  A = 
              x1         x2         x3         x4         x5
   x1      0.884    -0.1988     0.2346    -0.1524    -0.1915
   x2     0.1235     0.7245     0.3915    -0.1977    -0.1996
   x3    0.06063    -0.1976     0.6687     0.6427   -0.04188
   x4    0.03132    -0.0749    -0.1659     0.5566     -1.031
   x5   -0.00203   0.004721  -0.007366     0.2889     0.7387
 
  B = 
              u1
   x1   -0.05788
   x2    0.09314
   x3     0.0664
   x4     0.0368
   x5  -0.002879
 
  C = 
            x1       x2       x3       x4       x5
   y1  -0.3897  -0.4851   0.5155  -0.4198  -0.2795
 
  D = 
              u1
   y1  0.0001832
 
Sample time: 0.02 seconds
Discrete-time state-space model.

sysx = pade(sys3,nstates)

Warning: The model is already discrete. For discrete models, the "pade" command maps all delays to poles at z=0.

sysx =
 
  A = 
              x1         x2         x3         x4         x5
   x1      0.884    -0.1988     0.2346    -0.1524    -0.1915
   x2     0.1235     0.7245     0.3915    -0.1977    -0.1996
   x3    0.06063    -0.1976     0.6687     0.6427   -0.04188
   x4    0.03132    -0.0749    -0.1659     0.5566     -1.031
   x5   -0.00203   0.004721  -0.007366     0.2889     0.7387
 
  B = 
              u1
   x1   -0.05788
   x2    0.09314
   x3     0.0664
   x4     0.0368
   x5  -0.002879
 
  C = 
            x1       x2       x3       x4       x5
   y1  -0.3897  -0.4851   0.5155  -0.4198  -0.2795
 
  D = 
              u1
   y1  0.0001832
 
Sample time: 0.02 seconds
Discrete-time state-space model.

[DelayFreeSys,tau] = getDelayModel(sys3)

DelayFreeSys =
 
  A = 
              x1         x2         x3         x4         x5
   x1      0.884    -0.1988     0.2346    -0.1524    -0.1915
   x2     0.1235     0.7245     0.3915    -0.1977    -0.1996
   x3    0.06063    -0.1976     0.6687     0.6427   -0.04188
   x4    0.03132    -0.0749    -0.1659     0.5566     -1.031
   x5   -0.00203   0.004721  -0.007366     0.2889     0.7387
 
  B = 
              u1
   x1   -0.05788
   x2    0.09314
   x3     0.0664
   x4     0.0368
   x5  -0.002879
 
  C = 
            x1       x2       x3       x4       x5
   y1  -0.3897  -0.4851   0.5155  -0.4198  -0.2795
 
  D = 
              u1
   y1  0.0001832
 
Sample time: 0.02 seconds
Discrete-time state-space model.

tau =

  0x1 empty double column vector

tau

tau =

  0x1 empty double column vector

% Convert state-space model to transfer function
sys_tf  = tf(sys3);           
[n, d]  = tfdata(sys_tf, 'v'); 
isminphase([n,d])

ans = logical
   0

isallpass([n,d])

ans = logical
   0

% !Pre-filter, let's behave professionally, the original idea of such
% fitler is of Sam Chak from here and so I must cite him and credit his
% contribution: https://uk.mathworks.com/matlabcentral/profile/authors/5570092
%https://uk.mathworks.com/matlabcentral/answers/2112651-i-am-trying-to-find-the-lqr-controller-system-transfer-function-but-getting-error-as-arrays-have#answer_1449181
AdjGain = 0.1;
% with my additions of more scripting/coding control
fdst = length(d)

fdst = 6

fdend = length(d)

fdend = 6

fnst = length(n)-2

fnst = 4

fnend = length(n)

fnend = 6

Gf = AdjGain*tf(d(fdst:fdend),n(fnst:fnend),s_time)   

Gf =
 
               0.03949
  ----------------------------------
  0.02791 z^2 - 0.02041 z + 0.001796
 
Sample time: 0.02 seconds
Discrete-time transfer function.

% Display the transfer function of the Compensated System
disp('Transfer Function of the Compensated System:');

Transfer Function of the Compensated System:

tgain = length(A)

tgain = 5

%Desired controller:
res_agent= tgain*sys_tf

res_agent =
 
  0.0009161 z^5 - 0.0185 z^4 + 0.0812 z^3 - 0.1396 z^2 + 0.1021 z - 0.00898
  -------------------------------------------------------------------------
         z^5 - 3.573 z^4 + 5.559 z^3 - 4.663 z^2 + 2.088 z - 0.3949
 
Sample time: 0.02 seconds
Discrete-time transfer function.

figure
w = logspace(2,6,1000);
bode(res_agent,w),grid

figure
rlocus(res_agent)

Gcc = minreal(res_agent)        

Gcc =
 
  0.0009161 z^5 - 0.0185 z^4 + 0.0812 z^3 - 0.1396 z^2 + 0.1021 z - 0.00898
  -------------------------------------------------------------------------
         z^5 - 3.573 z^4 + 5.559 z^3 - 4.663 z^2 + 2.088 z - 0.3949
 
Sample time: 0.02 seconds
Discrete-time transfer function.

ms = minreal(Gcc);
zms = zero(ms)     % closed-loop zeros

zms = 
  14.8979 + 0.0000i
   2.6497 + 0.0000i
   1.2730 + 0.9131i
   1.2730 - 0.9131i
   0.1012 + 0.0000i

pms = pole(ms)     % closed-loop poles

pms = 
   0.6509 + 0.6264i
   0.6509 - 0.6264i
   0.7676 + 0.2623i
   0.7676 - 0.2623i
   0.7355 + 0.0000i

%ToDo: your conversion to angles:

%ToDo: add disturbances/noise simulation:

%Step response of the resultant system to step in amplitude 24.5 applied at
%time = 1 sec - 0.1 sec

%Evaluation:

tsim = 10;
setpointamp = 24.5; %Bar or Newton or mm
setpointapptime = 0.9; % sec
defaultsetpointpos = 0; % Bar or Newton or mm
Conf = RespConfig(Bias=defaultsetpointpos,Amplitude=setpointamp,Delay=setpointapptime)

Conf = 
  RespConfig with properties:

                Bias: 0
           Amplitude: 24.5000
               Delay: 0.9000
        InitialState: []
    InitialParameter: []

figure
step(Gcc,[0 tsim], Conf)
title('Plant+ FilteredLQR ™®© Controller Closed Loop step response')

isminphase(tfdata(tf(Gcc),'v'))

ans = logical
   0

isallpass(tfdata(tf(Gcc),'v'))

ans = logical
   0

stepinfo(Gcc)

ans = struct with fields:
         RiseTime: 0.0800
    TransientTime: 0.5200
     SettlingTime: 0.5200
      SettlingMin: 0.9336
      SettlingMax: 1.0794
        Overshoot: 5.1698
       Undershoot: 1.3945
             Peak: 1.0794
         PeakTime: 0.2600

figure
w = logspace(2,6,1000);
bode(Gcc,w),grid

figure
rlocus(Gcc)

figure
nyquist(Gcc)

figure
nichols(Gcc)

hasdelay(Gcc)

ans = logical
   0

td = totaldelay(Gcc)

td = 0

sysnd = absorbDelay(Gcc)

sysnd =
 
  0.0009161 z^5 - 0.0185 z^4 + 0.0812 z^3 - 0.1396 z^2 + 0.1021 z - 0.00898
  -------------------------------------------------------------------------
         z^5 - 3.573 z^4 + 5.559 z^3 - 4.663 z^2 + 2.088 z - 0.3949
 
Sample time: 0.02 seconds
Discrete-time transfer function.

sysx = pade(sys1,nstates)

Warning: The model is already discrete. For discrete models, the "pade" command maps all delays to poles at z=0.

sysx =
 
  A = 
               x1          x2          x3          x4          x5
   x1      0.9309     -0.3082     0.07682     -0.2736      0.1415
   x2     0.04795      0.9004      0.6453   -0.002722     -0.7354
   x3    0.006802    -0.07222      0.8497      0.7817     -0.4239
   x4    0.001494   -0.005394    -0.06559      0.6336      -1.242
   x5    0.000304  -0.0007171    -0.01521      0.2829      0.7553
 
  B = 
              u1
   x1   -0.05788
   x2    0.09314
   x3     0.0664
   x4     0.0368
   x5  -0.002879
 
  C = 
            x1       x2       x3       x4       x5
   y1  -0.3897  -0.4851   0.5155  -0.4198  -0.2795
 
  D = 
              u1
   y1  0.0001832
 
Sample time: 0.02 seconds
Discrete-time state-space model.

%[DelayFreeSys,tau] = getDelayModel(sys1)

%tau

%Constructed by
    %https://independent.academia.edu/PMazniker
    %[phone], https://join.skype.com/invite/oXnJhbgys7oW ™®©
    %https://diag.net/u/u6r3ondjie0w0l8138bafm095b
    %https://github.com/goodengineer
    %https://orcid.org/0000-0001-8184-8166
    %https://willwork781147312.wordpress.com/portfolio/cp/
    %https://www.youtube.com/channel/UCC__7jMOAHak0MVkUFtmO-w
    %https://nanohub.org/members/130066
    %https://pangian.com/user/hiretoserve/
    %https://substack.com/profile/191772642-paul-m

see my answer updated.

This is my original answer

I am on it, doing,

Better system identification will of course improve the controller(agent) design and implementation. LQ than can be integrated to MPC.

From what you uploaded the length of u and y is 135.

Can you make more measurement to catch all behaviour of your system, u from 0.3 to 23.5 with step with resolution 0.1 is so far is not enough. To tackle quantization noise also. Make more resolutions and span of value, let's 0.01 resolution and from 0.01 to 30.

Also the technique of system identification can be improved, let me get and I migh aid help making:

filter list squares based estimation or sys = ssest(u,y,nx,'Ts',Ts)

a few options I putted for you to choose most emended(optimal) for your systems.

Estimate ARX model using QR factorization.

Estimate state-space model using combined subspace method: MOESP algorithm for finding the matrices A and C, and N4SID algorithm for finding the matrices B and D.

Estimate state-space model using MOESP algorithm.

Estimate state-space model using N4SID algorithm.

Also the LQ I might do on paper I have more schemes and equations.

There are [g, x, l] = lqr (a, b, q, r, s, e) with more data for better lqr.

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

OK. Let's make it State of the Art, well-ordered to serve present and future real world natural goals, to control, the novel scientific heuristic I applied here is to advance the integrator (with mostly just found appropriate script built in expression of Sam Chak and integrated it with the LQR), make it more fitted to a specific signals and systems for industrial automatics, intelligent transportation, CNC machines more benign for humanity than aircrafts,

so q. N1: when it comes to implement the controller(agent) on real digital computer board processor on embedded C/C++ platform, the controller on mcu will read the actual plant(process, environment) y from sensor and A/D, while inside the C/C++ it will compute the controls(actions) u to go through transducer to the real physical plant(environment,process) and to the A,B,C,D, system or just to the real world plant? from K gains and A,B,C,D, matrices and precomputed already in Simulink, right?

While what with the needed for the controller(agent) state X, should it be computed inside the board in C/C++ during the controls(actions) generation in closed loop running mode from the precomputed and stored A,B,C,D, matrices and which? or must it (X) be estimated from the real world physical plant,- if so how and which adjustments to the control scheme can already help in real physical more efficient realizations and runs of it?

What if in reality the operating region is large and disturbances are harsh, such that the linearization and time variance corrections (due to tear, wear and exogenous factors) will be required each Nl steps?

What is the best directing stimulating material, visual datum for human brain to incite/induce/ most precise correctly grasping/apprehending and constructive, not destructive thinking process and full perception of all forms of signal at all pathes at only necessary points in time, re-transformations, effects on states during the regulation/control and setpoint commands towards full understanding of the control, regulation mathematics and closed loop plant-controller and input, output evolution operation flow subtleties?

Please, I need to receive, to get for upgrade of my RAM, Ethernet link, clockspeed and TCE NCE MPPL software, tools.

[Mathieu NOE]

@Wiktoria

hello again - so I looked at the identification data you shared recently

you did not mentionned at which sample rate that was acquired , but as you said in another comment :

The actuator itself doesn't have such a big delay at all! Actually, the delay is around 20 ms only.

I assumed that the sample rate for this step (identification) acquisition is 1 ms

Am I right ?

is there any confusion between a delay in samples (as we can clearly see 20 samples) or in ms ?

then there is probably a correction to be done when you identified a SS model that is supposed to be valide with a 20 ms sampling rate (vs 1 ms in the identification data collection)

[Wiktoria]

@Mathieu NOE Actually this identification was acquired with 0.01 sampling rate. I have data from external data logger and when downloading data from it, I can choose what sampling rate I want. I have previously identified system with 0.02 sampling, but then tried experimenting. I can also add data with 0.02s sampling rate, but when I tested some models I had, the previous ones were unsatisfactory.

When I mentioned the actuator, I heard before that a servomechanism's delay can be around 20 ms and it is a mechanical and electrical delay. That delay was visible in logs too - the reaction to impulse wasn't instant, servo needed time to move, it was fast, but around 10 to 20 ms depending on the range of the move.

Now that you mentioned it, I am confused with the identification. It was my first time doing it. I'm unsure how it should be done. In system identification toolbox, the only thing I could put was sampling rate of my data, input and output.

[Wiktoria]

@Pavl M. My sampling rate isn't 0.1s, it is 0.01s. Currently it appears to me that the model itself is slow, if it even is possible. I have changed control and gains, my system appears to be much faster in simulation but on plant it still has a delay, although smaller. I will check the solutions you offered with the plant and wll give you feedback if LQG was better. Thank you!

[Mathieu NOE]

@Wiktoria

ok I let you check what was the actual sampling rate during identification

when you exported you data in the mat file , did you make any rate conversion between the data logger and mat file ?

to sum up , the 200 ms delay that appears in your SS model seem indeed incorrect vs reality

[Wiktoria]

@Mathieu NOE The data logger resamples data according to the user's input - when downloading data, I have to input the sampling I want. If I don't do that, the sampling is set up automatically and it's usually 0.01s. I have two sets of data prepared - one is sampling every 0.01s and the other is 0.02s. I will share the data with 0.02s sampling in the question above.

I also wonder if the data I'm using is correct in the first place, meaning that the impulse that I gave didn't consist of 1 bar - I have set a requested brake pressure to 24.5 bars and the step went from 0 to said pressure, as it's visible in data I have shared. What makes me wonder is that I have heard that when identifying a plant, the most common thing to do is to use a unit step. Is there any need for rescaling, or when using System Identification Toolbox, MATLAB does that by itself already? Or perhaps it's not necessary at all?

My other question is: if my model has such a delay, I believe I should make another one. How can I avoid making such a mistake again? Meaning that I was mostly just following one of the guidelines when it comes to identification and I don't know where such a delay was added in the process. Perhaps it's an issue with how data is prepared for identification?

[Mathieu NOE]

@Wiktoria

the goal here is to identify the "plant" that is the portion of the system that goes from your control board output to its input (where the pressure sensor is connected)

so you send a unitary step "software" signal that is then converted to a physical voltage (this is done by the DAC of your control board, or it can also be a PWM signal) that goes into the valve drive , which itself control the pressure in your hydraulic cylinder , ...etc ... etc

I believe I am not telling you here anything you are not already aware, but this is just to put hings in place again (helps me too to construct my answer)

so when we say we generate a unitary step , this is only inside your control board , at the software level. Your DAC may have a certain gain (how much volts at the output connector for a unitary step ?) , and same for the other devices in your system, each one has a different DC gain (and poles /zeroes).

so if I were you and someone gives me that "plant" to identify and I have no clue about what's in there , I would go step by step :

• what kind of test signal is appropriate for this task . Not all systems react linearly with steps or impulses. The physics of the "DUT" (device under test) should help decide which signal is best for the job. We want a "plant friendly" test signal , so that we get the desired amount of data to perform a good identification, while remaining in the physical limits and if possible in the linear domain. If you go too fast , your system may react with a certain slope due to its physical velocity limit (this is especially true with hydraulic actuators). If your test signal amplitude is too low , you get mostly noise , if too high you may hit the end stops of reach velocity limit
• so when confronted with a new system I like to start very gently (low amplitude signals) and then adapt the test signal characteristics according to how the output looks like. So for the system you have described , step is good IMO, but the need to do a full swing (max amplitude) input may be questionable. I would have probably done several attemps with increasing amplitudes and check whether the output is proportionnal to input or if we start to see some non linearity (staturation) . I probably start with looking at the signals on the scope at diferent amplitudes (and from there I should already know what are the typical characteristics of my system) , and once I have narrowed down to what I really want then I would do the official records that I will use for identification purpose.
• it may look a bit overkill and waste of time , but when you discover a new system it may actually save you time , because the more you know about the system and it's limits , the better you can refine your identification process and avoid going forth and back because you need to re-check your data.

if a unitary step generates a 24 bars response (and the system remains in the linear operating range) , there is no difference in term of system identification result vs having a smaller step amplitude calibrated for a 1 bar response (but I assume the shape of the output may be slightly different even rescaled). I am not sure which matlab function you used , but you may have to rescale the signal you feed to the function if you are using a non unitary step signal (you probably have to revisit the documentation and make some tests with a simple case model - or demo code)

also , it would be nice if in your record we could see when the step signal switch from zero to one (because today I only see the pressure signal after the transition, so I may have a doubt when the transition actually occured).

sorry for saying something you know already , both signals must be recorded synchronously

anyway , if you could share some records, I'll be glad to have a look at them and give you some feedback (no pun intended !)

all the best !

[Wiktoria]

@Mathieu NOE Hello,

I have just shared the data with the transition, sampling 0.02s, just like before.

Regarding what you have said, I will try different amplitudes (just like you adviced) to try to figure the system out. The main issue though is that there is a lot of noise on the brake sensors and making a unitary step (so for going from 0 to 1) shows almost no change. I will try with smaller steps though and will see if the response changed, and I'll try to identify a plant yet again.

Also, today I tested the controller on my plant and it is way faster than before, but it still has a certain delay - probably the same one that is also visible in identified model. Do you think there is any way to fix it? Perhaps a different state space equations or Smith predictor? Additionally, there appears to be a much larger delay in the very first moment of regulation than in all the others. It appears to me that after I initialize the controller, the request takes a longer time to process (at least that's my guess) and a much larger delay appears. All the responses later on, after the first move, are much faster. The request is sent via CANbus, there shouldn't be much delay there. I wonder if that's also an issue of my system. I have shared the response of my controller too (screenshots in pdf file), so if you would like to look at that too, I would be very thankful.

Thank you very much for your help and hope you are having a wonderful day!

[Mathieu NOE]

hello @Wiktoria

both new data files seems to show a delay around 7 samples ( x 0.02s = 140 ms, so not very far from your initial statement about 200 ms delay)

with_sampling_0.02.mat :

transition.mat

what I see on the graphs in your pdf document is a bit surprising to me as the system (in closed loop) seems to show longer and more variable delays when transition from low to high pressure input - a contrario when going downwards from high to low pressure input , it's much faster. is there somewhere a non symmetry in your plant that explains why it's reacting much faster in one direction vs the other ?

here below we can see the first transition (low to high) takes almost 1 s , whereas the next ones are much faster ...

[Wiktoria]

@Mathieu NOE Well, there shouldn't be any non-symmetrical behaviors, I actually have no idea why there is such a large delay when servo is letting go of the pedal. I wonder if there could be any hidden non linear behavior in the plant or in the system?

I am going to perform the identification again and see if any of this behavior gets changed. My current biggest guess is that either that's the identification fault or somehow the plant is non-linear in the lower pressure values.

Answer 4

[Pavl M.]

hmm... quite interesting Interesting... especially if you have the delay large only on first activation and then the delay is reducing... alpha it is needed to understand fully exactly why and what is the only correct sources of the delay, is delay caused before the pedal pulling end or after it?

let me think, see... there are many work to do on this project ... firstly it must needs that you will be able to control braking robust in harsh real world automotive environments (take into accound disturbances there) via CAN bus or any other digital bus (b.t.w. there are alternatives to CAN with the same goal and maybe better, it is needed to research as well) to your mcu board and then from board through servo, pedal and all the pneumohydraulic braking hoses till the real wheels ... ( next real work scope for this project I see a) to make more algebra, eigen calculus and matrix math to further optimize the control(regulation) strategem (let me do it, how much it costs among the rest?), b) reveiw the automotive braking system hardware designs practices, to find most suitable for your(also the relative cost of this part?), c) understand more system identification algorithms as much deliberately and not by trial and error, agreable only just estimate how much it takes to read more than 10 conventional book pages per each function in documentation, large documentation for this is not undreasonable, (there are more than 4 so far I found, also from frequency domain data system can be estimated), understand non linearities and time variations there, understand how to estimate general parametric and time varying non-linear system as in nature, understand more controlers (regulators, agents, compensators) succesfull, real world designs and realistic p.o. that we do.

From recent system identification data slicot wrapped TCE NCE MPPL estimated now the system better and without many previous rank defficiency warnings. The cost of control is not unreasonably high. There are no free lunches of course. These are results attached for your recent datum for system input, output passed through 4 system indentification algorithms with 0.01 and 0.02 sampling rates and estimated for 7, 5 and 3 states so far with Kalman and covariance matrices. I have for 0.02 your samples also 7 states estimated and for 0.01 5 states in .zip archives, let me know whether to send ye. For real-time applications, you can generate code from the Multistage Nonlinear MPC block.

The system must work under hard real time constraint, so the less delay the better and so the controller must not be very complex, the simplest one is just to find and cancel zero (numerator parts accordingly) in System Identified transfer function which causes the delay. Second it must be prone/robust to noise and disturbances, which increases control energy demand by some margin.

Why just 107 samples, why the samples and not others, can you make more samples, 300?

A. I understand your concern about the delays and controller other than LQR, like Smith Predictor to fasten/rapidize the system and remove dead times, but whether the observed delay is a dead time, I think no. Kindly correct me if m wrong. I have worked today on it and there are left to complete a few from my side. See attached script to this answer. TCE NCE MPPL Matlab detected no inherent delays in the systems.

Mathematically, delay in frequency domain / transfer function for continuous time case is *exp(-tau*s), which is itself a long polynom (see Maclauren series = 1 - tau*s + tau*s)^2 + ...) or *z^(-tau) for discrete time case, which is also can be approximated as Taylor series in z.

so to assume the delay must be already resembled after System Identification (after Combined Subspace Methode, SSSID or Kalman like estimator may be employed in System Identification) process in the most valued precise A,B,C,D, matrices.

and I opened your with_transition.dat and with_sampling_002.dat and found that the less delay ( both are approximatelly 20 ms) is in with_transition.dat, the longes delay in you initial identification_data.dat .

Let me know do yoy need it in Julia, Pyton, R (high level design), C/C++ or Python for embedded for low level design?

Also can we do the embeddable to MCU board firmware right from script (novel RT embedded software and firmware development technique)? Is it better(and accodring to which criteria of the engineering design and contrustion technique evaluation?) to do firmware from script commands or from Simulink like block scheme / diagram?

My answer the Full State Feedback LQR Controller with pre-filter (select which one is better for you (r-K*X)*filter or (r*filter - K*X) ) designed and applied(simulated) shown significant reduction in delay (rise time from 0.18 in open loop system to 0.08 sec in closed loop, settling and transient time from 1.16 to 0.52).

Of course it is yet not perfect, specifically not very robust, non minimal phase, multivariate can be added to it.

Should I put more care for the delays or is the Full State Feedback LQR Controller with pre-filter is sufficient enough to reduce your delays so far? Please let me know. There is also MataveID contains realization identification, polynomal algorithms and subspace algorithms, which needs practical examples and the codes are at github, if they will be willing to pay salary in exchange for , w/o further ado, no strings attached, can do it, while here who will be needing not open source software?

Which type of delay is it input, output or internal delay? Mathematically, whether the delay is fully expressed by the System Identification toolbox provided transfer function zeroes (nominator coefficients)?

Next matlab functions are useful for some in-depth analysis (sys1 is the discrete time system got from System Identification process), why they detect no internal or input/output delays in your system, they display yet no input nor output no internal dealy there, what is the gist of these:

stepinfo(sys1)

hasdelay(sys1)
td = totaldelay(sys1)
sysnd = absorbDelay(sys1)
sysx = pade(sys1,nstates)

[DelayFreeSys,tau] = getDelayModel(sys1)