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

115 次查看（过去 30 天）
Wiktoria
Wiktoria 2024-11-14，20:41
编辑： Pavl M. 3 minutes 前
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;
  8 个评论
显示 6更早的评论
Pavl M.
Pavl M. 2024-11-17，17:32
编辑：John Kelly 2024-11-18，15:31
Attn.:
The original idea to employ the Hermittian 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
X_dot = A*X+B*U,
then the control output can be computed as:
LQR:
U = (setpoint - y_from_real_phys_system)*Klqr(N+1) - (Klqr1.. KlqrN)*(X1 .. XN) (designed so far with augmented)
U = (setpoint - y_from_real_phys_system)*Klqr(1) - (Klqr2.. KlqrN)*(X2 .. XN)(not augmented)
or U = F1(setpoint - y_from_real_phys_system) - (Klqr1...KlqrN)*(X1..XN)
or (see chk1.m script uploaded)
LQI:
u = DiscreteIntegratorApproximator(r - y_from_real_phys_system)*Klqi(N+1) - Klqi(1..N)*X(1..N) (most precise)
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)
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（3 个）

Mathieu NOE
Mathieu NOE 2024-11-15，11:51
Ran in:
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
  2 个评论
Mathieu NOE
Mathieu NOE 2024-11-15，12:03
Ran in:
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.
Pavl M. 2024-11-15，18:35
编辑：Pavl M. 2024-11-15，18:37
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.

请先登录，再进行评论。

Pavl M.
Pavl M. 2024-11-15，11:06
编辑：Pavl M. about 14 hours 前
It is indeed really factually very valued, yet undervalued and quite intersting, utile question.
(Updated, 17 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)
home
clear all
close all
pack
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:')
isstable(sys1)
isminphase(tfdata(tf(sys1),'v'))
isproper(sys1)
Qc= ctrb(A,B)
controllab = rank(Qc)
x = 1;
disp('The pair (A,B) must be stabilizable')
ic = [0 0 0 0 0];
figure
step(sys1);
Gd1 = c2d(sys2,s_time,'impulse')
Gd2 = d2d(sys1,s_time,'tustin')
%figure
%step(Gd1)
figure
step(Gd2)
[A_size_row, A_size_column] = size(A);
Q = 70*eye(A_size_row+1)
R = 7 %*eye(ninputs);
R2 = 5 %*eye(ninputs+1);
N = ones(nstates,1)
N2 = ones(nstates+1,1)
%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
%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)
%% For Just LQR:
[K1,S1,P1] = lqr(sys1,Q2,R)
%% For just LQI:
[Klqi,Slqi,elqi] = lqi(sys1,Q,R2,N2)
%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)
%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
H2 = A - B*K2
%% 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);
% Convert state-space model to transfer function
sys_tf = tf(sys3);
[n, d] = tfdata(sys_tf, 'v');
% Pre-filter
Gf = -0.1*tf(d(end), n,s_time)
% Display the transfer function of the Command Compensated System
disp('Transfer Function of the Command Compensated System:');
Gcc = minreal(Gf*sys_tf)
ms = minreal(Gcc);
zms = zero(ms) % closed-loop zeros
pms = pole(ms) % closed-loop poles
%Step response of the resultant system to step in amplitude 12 applied at
%time = 1 sec
% replaces loop:
tsim = 10;
setpointamp = 24; %Bar or Newton or mm
setpointapptime = 1; % sec
defaultsetpointpos = 0; % Bar or Newton or mm
Conf = RespConfig(Bias=defaultsetpointpos,Amplitude=setpointamp,Delay=setpointapptime)
figure
step(Gcc,[0 tsim], Conf)
title('Plant+Controller LQR Closed Loop system step response')
%% PID
P1 = tf(sys1)
[C, info] = pidtune(P1, 'PIDF2')
C2tf = tf(C)
Cr = C2tf(1)
Cy = C2tf(2)
Gclp = Cr*feedback(P1,Cy,+1)
Gcl = minreal(Gclp);
zcl = zero(Gcl) % closed-loop zeros
pcl = pole(Gcl) % closed-loop poles
% 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
%[K,S,e] = lqry(sys,Q,R,N)
%% More robust gaussian servo controller:
nx = 5; %Number of states
ny = 5; %Number of outputs
%Disturbance parameters, covariance matrix, here to tune the QXU,QWV,QI (find most
% optimal values):
Qn = Q2+100*[0.04 0.02 0 0 0; 0.02 0.01 0 0 0; 0 0 0.01 0 0; 0 0 0 0 0; 0 0 0 0 0]
Rn = 0.1;
R = sqrt(Qn); %[0.01 0 0 0 0; 0 0.02 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0]
QXU = blkdiag(eye(nx),R)
QWV = blkdiag(Qn,Rn)
QI = 1 ; %eye(ny)
%Servo-controller design (Attn.: EXPERIMENTAL,IT'S LEFT TO TUNE THE MATRICES):
KLQG1 = lqg(sys1,QXU(1:A_size_row+1,1:A_size_column+1),QWV,QI,'1dof')
KLQG2 = lqg(sys1,QXU(1:A_size_row+1,1:A_size_column+1),QWV,QI,'2dof')
%Constructing closed-loop system after designing:
%Ther are many similar questions this answer will fit(answer) to:
%https://www.mathworks.com/matlabcentral/answers/484846-how-to-use-the-generated-lqg-controller
C=KLQG1;
P = sys1;
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)
figure
step(clsys_m,[0 tsim], Conf)
title('KLQG1 closed loop system response')
C=KLQG2;
P = sys1;
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)
figure
step(clsys_m2,[0 tsim], Conf)
title('KLQG2 compensated closed loop system response')
%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
%+380990535261, 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.
  2 个评论
Wiktoria
Wiktoria 2024-11-15，13:01
编辑：Wiktoria 2024-11-15，13:14
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.
Pavl M. 2024-11-16，7:11
编辑：Pavl M. 2024-11-17，19:34
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_dot = A*X+B*U,
then the control output can be computed as:
LQR:
U = (setpoint - y_from_real_phys_system)*Klqr(N+1) - (Klqr1.. KlqrN)*(X1 .. XN) (designed so far with augmented)
U = (setpoint - y_from_real_phys_system)*Klqr(1) - (Klqr2.. KlqrN)*(X2 .. XN)(not augmented)
or U = F1(setpoint - y_from_real_phys_system) - (Klqr1...KlqrN)*(X1..XN)
LQI:
u = DiscreteIntegratorApproximator(r - y_from_real_phys_system)*Klqi(N+1) - Klqi(1..N)*X(1..N) (most precise)
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)
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]?
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))?

请先登录，再进行评论。

Paul
Paul 2024-11-17，17:14
Ran in:
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.
  5 个评论
显示 3更早的评论
Mathieu NOE
Mathieu NOE 2024-11-18，9:49
@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.
Pavl M. 2024-11-18，18:03
编辑：Pavl M. 3 minutes 前
@Wiktoria,
From where they took next quadratic cost:
H = a*y^2 + b*yi^2
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 reduced dealy by reducing t_rise(rise time) and settling_time is already not bad, how phase and gain margins performance characteristics resemble the latency reduction/quantity?
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:
nx = 5; %Number of states
ny = 5; %Number of outputs
%Disturbance parameters, covariance matrix, quite small as in real world:
Qn = Q2+100*[0.04 0.02 0 0 0; 0.02 0.01 0 0 0; 0 0 0.01 0 0; 0 0 0 0 0; 0 0 0 0 0]
Rn = 0.1;
R = sqrt(Qn); %[0.01 0 0 0 0; 0 0.02 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0]
QXU = blkdiag(eye(nx),R)
QWV = blkdiag(Qn,Rn)
QI = 1 ; %eye(ny)
%Servo-controller design:
KLQG1 = lqg(sys1,QXU(1:A_size_row+1,1:A_size_column+1),QWV,QI,'1dof')
KLQG2 = lqg(sys1,QXU(1:A_size_row+1,1:A_size_column+1),QWV,QI,'2dof')
%Constructing closed-loop system after designing:
%Ther are many similar questions this answer will fit(answer) to:
%https://www.mathworks.com/matlabcentral/answers/484846-how-to-use-the-generated-lqg-controller
C=KLQG1;
P = sys1;
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)
figure
step(clsys_m,[0 tsim], Conf)
title('KLQG1 closed loop system response')
C=KLQG2;
P = sys1;
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)
figure
step(clsys_m2,[0 tsim], Conf)
title('KLQG2 compensated closed loop system response')
(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
sys = ssest(data,nx) estimates a continuous-time state-space model using the time-domain or frequency-domain data in the data object data.
Then you can get more data for your A,B,C,D matrices, ushc us nx = 5 or 4 or more.
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.
ssIt 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?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 State-Space Control Design 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by