Scribd
Upload
42 views

Check Stability of A: Appendix MATLAB Code

The MATLAB code provides 4 sections that: 1. Check the stability of system matrix A and finds it is unstable 2. Checks the controllability of system AB and finds it is controllable 3. Uses LMI toolbox to find optimal state feedback gain K to stabilize the system 4. Plots the closed loop system states over time using the optimal controller gain K

Uploaded by

Apra Saurabh Gupta
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
42 views

Check Stability of A: Appendix MATLAB Code

The MATLAB code provides 4 sections that: 1. Check the stability of system matrix A and finds it is unstable 2. Checks the controllability of system AB and finds it is controllable 3. Uses LMI toolbox to find optimal state feedback gain K to stabilize the system 4. Plots the closed loop system states over time using the optimal controller gain K

Uploaded by

Apra Saurabh Gupta
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
You are on page 1/ 6

APPENDIX

MATLAB Code:

1. Check stability of A

clear all;
clc;

% system description
A = [0 1 0 0; 20.601 0 0 0; 0 0 0 1; -0.4905 0 0 0];
disp(eig(A)); % Calulate and display the Eigen Values of A

Solution:
A= 0 1.0000 0 0
20.6010 0 0 0
0 0 0 1.0000
-0.4905 0 0 0

eig(A) = 0
0
4.5388
-4.5388

Matrix A is Unstable

2. Check Controllability of AB

clear all;
clc;

% system description
A = [0 1 0 0; 20.601 0 0 0; 0 0 0 1; -0.4905 0 0 0];
B = [0; -1; 0; 0.5];

Co = ctrb(A,B); %Compute controllability matrix

unco = length(A) - rank(Co) % Determine the number of uncontrollable


states

Solution:
unco = 0

Number of Uncontrollable states is NIL hence AB is Controllable


3. Optimal Solution using LMI toolbox

clear all;
clc;

% system description
A = [0 1 0 0; 20.601 0 0 0; 0 0 0 1; -0.4905 0 0 0];
B = [0; -1; 0; 0.5];

% performance matrices
Q = [200 0 0 0; 0 2 0 0; 0 0 2 0; 0 0 0 2];
R = 2;

% LMI flow
setlmis([])
X = lmivar(1, [4,1]); % Decision variable

lmiterm([-1 1 1 X],1 ,1); % X>0


lmiterm([2 1 1 X],A ,1, 's'); % AX +XA'
lmiterm([2 1 1 0],-B*inv(R)*B'); % -B*inv(R)*B'
lmiterm([2 1 2 X],1 ,1); %X
lmiterm([2 2 2 0],-inv(Q)); % -inv(Q)

LMIs = getlmis;
[TMIN,XFEAS] = feasp (LMIs); % feasibility of LMI (TMIN<0 implies feasible)

X = dec2mat (LMIs, XFEAS, X); % numerical form of decision variable


P=inv(X); % transforming back to P

K = inv(R)*B'*P % state feedback gain

D = (A-B*K); % closed loop matrix


eig(D) % closed loop eigenvalues

Solution:
Solver for LMI feasibility problems L(x) < R(x)
This solver minimizes t subject to L(x) < R(x) + t*I
The best value of t should be negative for feasibility
Iteration : Best value of t so far

1 0.030747
2 0.015314
3 0.015314
4 3.103825e-03
5 3.103825e-03
*** new lower bound: -0.080638
6 3.103825e-03
*** new lower bound: -0.014184
7 6.389709e-04
8 -5.522330e-04
*** new lower bound: -0.011082

Result: best value of t: -5.522330e-04


f-radius saturation: 0.000% of R = 1.00e+09

K = -80.1449 -17.7407 -2.3320 -6.2538

Eig (A-B*K) = -8.6488 + 0.0000i


-4.6197 + 0.0000i
-0.6727 + 0.3465i
-0.6727 - 0.3465i

System has stable closed loop poles, hence K is optimal Controller gain

4. Plot of States with respect to time

Clc;
clear;
close all;

% System matrix
A = [0 1 0 0; 20.601 0 0 0; 0 0 0 1; -0.4905 0 0 0];

% Control matrix
B = [0; -1; 0; 0.5];

% Controller gain
K = [-80.1449 -17.7407 -2.3320 -6.2538];

% Output matrix
C = [1 0 0 0;0 1 0 0; 0 0 1 0; 0 0 0 1];

% Feed forward matrix


D = [0;0;0;0];

% Initial state
x0 = [0.1; 0; 0; 0];

% State space modelling with controller feedback


sys_cl = ss(A-B*K,B,C,D);

t = 0:0.01:15; % Defining Time axis


[y,t,x] = initial (sys_cl,x0,t);

figure
subplot(2,2,1) % Plot of Angular Position of Pendulum w.r.t.time
plot(t,x(:,1),'k') % Plot color to be black
xlabel('Time (t)')
ylabel('\theta(t)')
title('Angular position')
grid on
set(gca,'FontSize',12)

subplot(2,2,3) % Plot of Angular Velocity of Pendulum w.r.t.time


plot(t,x(:,2),'b') % Plot color to be blue
xlabel('Time(t)')
ylabel('$\dot{\theta}(t)$', 'Interpreter','latex')
title('Angular Velocity')
grid on
set(gca,'FontSize',12)

subplot(2,2,2) % Plot of Position of Cart w.r.t.time


plot(t,x(:,3),'m') % Plot color to be magenta
xlabel('Time (t)')
ylabel('x(t)')
title('Cart poition')
grid on
set(gca,'FontSize',12)

subplot(2,2,4) % Plot of Velocity of Cart w.r.t.time


plot(t,x(:,4),'g') % Plot color to be green
xlabel('Time (t)')
ylabel('$\dot{x}(t)$', 'Interpreter','latex')
title('Cart velocity')
grid on
set(gca,'FontSize',12)

figure
plot(t,x) % Plot of all states w.r.t.time
xlabel('Time')
ylabel('States')
legend('x_1','x_2','x_3','x_4')
set(gca,'FontSize',18)

You might also like