How can I reduce oscilltions in control signal in LQI controller and why response is so similar even using different models?
Show older comments
I am trying to control a servomechanism that creates pressure by pulling a brake pedal. I have identified the system using System Identification Toolbox. I got the state space matrices, implemented an LQI controller and although the system works, the control signal oscillates when the impulse comes, even when changing reference signal from for example 5 bars to 0.
Although the oscillations in the beginning are small, I want to reduce them as much as possivle, since if the reference signal gets bigger, so do these oscillations. I have tried changing the model I am using and I tried to identify the system again, but surprisingly, this seems to always happen. I am confused why, when I change a model, the behavior is the same. Of course I tried changing the weight matrices Q and R, but even when the oscillations got lower, time of regulation got significantly longer and I want to avoid that too. I want to implement it "step by step", as I want to rewrite it later on to C.
I am completely new to control theory and I wonder what I am doing wrong, so I will be grateful for any help. My current code:
A = [-0.4628, 0.3217, 0.5789, -0.3823, -0.0507;
0.7593, 1.0396, -0.2008, -0.2438, -0.0678;
0.2106, 0.6490, 0.3973, 0.8735, 0.2070;
0.3128, 0.3380, 0.1524, 0.0914, -0.1088;
0.0353, 0.0419, -0.0355, 0.0682, 0.7576];
B = [0.0026;
0.0088;
-0.0036;
0.0490;
-0.1443];
C = [1.1099, 1.3983, 0.8999, -0.0136, 0.0652];
D = 0;
sys = ss(A, B, C, D, 0.01);
step(sys);
[A_size_row, A_size_column] = size(A);
Q = 4*eye(A_size_row+1);
R = 1;
A_aug = [A, zeros(A_size_column, 1); -C, 0];
B_aug = [B; 0];
[K, ~, ~] = lqr(A_aug, B_aug, Q, R);
K_i = K(6);
dt = 0.02;
T = 100; %in miliseconds
time = 0:dt:T;
x = [0;
0;
0;
0;
0];
x_i = 0;
y = [];
u = [];
ref = [];
states = [];
r=0;
for t = time
if t == 50
r = 5;
end
u_t = -K_x * x - K_i * x_i;
x_dot = A * x + B * u_t;
x = x + x_dot * dt;
y_t = C * x;
x_i_dot = r - y_t;
x_i = x_i + x_i_dot * dt;
y = [y, y_t];
u = [u, u_t];
ref = [ref, r];
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;
1 Comment
Sometimes you cannot achieve everything you want via a particualar controller structure. So sometimes it is a good idea to play around with how you set your reference. So oneway to deal with it is, since your system is slow enough, you can apply a ramp input that reaches to the end value within 10 seconds instead of step:
A = [-0.4628, 0.3217, 0.5789, -0.3823, -0.0507;
0.7593, 1.0396, -0.2008, -0.2438, -0.0678;
0.2106, 0.6490, 0.3973, 0.8735, 0.2070;
0.3128, 0.3380, 0.1524, 0.0914, -0.1088;
0.0353, 0.0419, -0.0355, 0.0682, 0.7576];
B = [0.0026;
0.0088;
-0.0036;
0.0490;
-0.1443];
C = [1.1099, 1.3983, 0.8999, -0.0136, 0.0652];
D = 0;
sys = ss(A, B, C, D, 0.01);
% step(sys);
[A_size_row, A_size_column] = size(A);
Q = 4*eye(A_size_row+1);
R = 1;
A_aug = [A, zeros(A_size_column, 1); -C, 0];
B_aug = [B; 0];
[K, ~, ~] = lqr(A_aug, B_aug, Q, R);
K_i = K(6);
K_x = K(1:5);
dt = 0.02;
T = 100; %in miliseconds
time = 0:dt:T;
x = [0;
0;
0;
0;
0];
x_i = 0;
y = [];
u = [];
ref = [];
states = [];
r=0;
for t = time
if t >= 50 & t <= 60
r = (t-50)*.5;
elseif t>= 60
r = 5;
else
r = 0;
end
u_t = -K_x * x - K_i * x_i;
x_dot = A * x + B * u_t;
x = x + x_dot * dt;
y_t = C * x;
x_i_dot = r - y_t;
x_i = x_i + x_i_dot * dt;
y = [y, y_t];
u = [u, u_t];
ref = [ref, r];
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;
ylim([-10 10])
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;
Accepted Answer
Mr. Pavl M.
on 28 Oct 2024
Edited: Mr. Pavl M.
on 28 Oct 2024
Good day.
The topic, theme, the agenda of the special regulator control strategy design, verification and implementation for custom actuators and extended pilot control is no doubt very valuable, interesting and I commenced shaping, sharpening, further enhancing it within the process in progress.
The scope at hand actually opens quite valuable large field control area of betterment workouting.
Let's make optimal controller for it,
What are the constraints on control effort(u) energy and speed?
What are the required settling time, overshoot, rise time?
As soon as I complete, contact me for results.
So far solved:
Found:
Whether the system is stabilizable:
In MPPL NCE TCE:
disp('Whether the system is stabilizable:')
isstabilizable(sys)
isstable(sys)
disp('The pair (A,B) must be stabilizable')
ans = 1
ans = 1
The pair (A,B) must be stabilizable.
The system is stable and stabilizable and so optimal Ricatti control exist.
Just to find it by finding the optimal Ki and Kx
1. In your code you actually activated Linear Quadratic Regulator (LQR) lqr function and not Linear Quadratic Integral (LQI) lqi methode.
2. For what to augment the matrices?
3. Try also good: reg = lqg(sys,QXU,QWV,QI,'1dof')
4. Why not to use just PID or fractional PI controller for pressure control?
5. Let me make also Simulink model for it:
It for sure widens and limits the solution space to find the optimal controller for large familty of plants.
Upon working on your code, the osciallation alleviated, just large signals appeared to fix.
So far I've constructed Gaussian Servo 1DOF controller clsys. Hope this to continue to contribute positive value to solving.
Also in this code, which is just a begining of more valuable commercial project I propose to create from it, there is no oscillations in control signal:
I created Simulink structure for it. I think I can do it less oscillatory, more stable. It must needs to get necessary support investment from industry to me, to make it more. Listen, see, reasonable: For more elegant, precise and professional completed results (I can do if you hire / pay me) more substantial funds investment is necessary
clc
clear all
close all
format long
%pkg load control
A = [-0.4628, 0.3217, 0.5789, -0.3823, -0.0507;
0.7593, 1.0396, -0.2008, -0.2438, -0.0678;
0.2106, 0.6490, 0.3973, 0.8735, 0.2070;
0.3128, 0.3380, 0.1524, 0.0914, -0.1088;
0.0353, 0.0419, -0.0355, 0.0682, 0.7576];
B = [0.0026;
0.0088;
-0.0036;
0.0490;
-0.1443];
C = [1.1099, 1.3983, 0.8999, -0.0136, 0.0652];
D = 0;
dt = 0.02;
aug_level = 0;
gain = 2;
sys = ss(A, B, C, D, dt,'InputName','u','OutputName',{'y'});
figure
step(sys);
%isstabilizable(sys)
%isprop(sys)
isstable(sys)
disp('The pair (A,B) must be stabilizable')
[A_size_row, A_size_column] = size(A);
Q1 = gain*eye(A_size_row+aug_level+1);
Q2 = gain*eye(A_size_row+aug_level);
R = 1; %eye(size(A)) requires square of size 1;
if aug_level == 0
A_aug = A;
else
A_aug = [A, zeros(A_size_column, aug_level); -C, zeros(1,aug_level);zeros(aug_level-1,A_size_column+aug_level)]
end
B_aug = [B; zeros(aug_level,1)]
[K1,XG1,L1] = lqi(sys,Q1,R)
%For Continuous system:
%[X2, L2, G2] = care(A_aug,B_aug,Q2,R)
%for discrete system:
[X3,L3,G3] = dare(A_aug,B_aug,Q2,R)
[K,XG,L] = dlqr(A_aug, B_aug, Q2, R)
%Optional!, not final, ToDo: in progress
K_i = -min(X3)
K_x = -X3
T = 100; %in miliseconds
time = 0:dt:T;
x = [0;0;0;0;0];
x_i = x;
y = [];
u = [];
ref = [];
states = [];
r=0;
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;
x = x + x_dot * dt;
y_t = C * x;
x_i_dot = r - y_t;
x_i = x_i + x_i_dot * dt;
y = [y, y_t];
u = [u, u_t];
ref = [ref, r];
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;
%N = eye(size(A))
%[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, quite small as in real world:
Qn = [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.01;
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(0.1*eye(nx),R)
QWV = blkdiag(Qn,Rn)
QI = 1 ; %eye(ny)
%Servo-controller design:
KLQG1 = lqg(sys,QXU(1:A_size_row+1,1:A_size_column+1),QWV,QI,'1dof')
KLQG2 = lqg(sys,QXU(1:A_size_row+1,1:A_size_column+1),QWV,QI,'2dof')
%Plot outcome - control product response:
figure
step(KLQG1)
title('1DOF Servo Controller designed')
%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 = sys;
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)
title('Constructed closed loop system response')
x = [0;0;0;0;0];
x_i = x;
y = [];
u = [];
ref = [];
states = [];
r=0;
%Todo: simulate it in loop:
%noise = 0.001;
%for t = time
% if t == 50
% r = 5;
% end
% u_t = ? + noise;
% x_dot = clsys.A * x + clsys.B * u_t + noise;
% x = x + x_dot * dt;
% y_t = KLQG1.C * x + KLQG1.D*u_;
% x_i_dot = r - y_t;
% x_i = x_i + x_i_dot * dt;
% y = [y, y_t];
% u = [u, u_t];
% ref = [ref, r];
% 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(sys,Qn,Rn,Nn);
%Connect the and sys to form the working loop, it should bring more performant control:
%K4 = lqi(sys,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')
More Answers (0)
Categories
Find more on Linear Model Identification in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
