SIMULINK CONTROLLER DESIGN

Goal 1: Form a-bk from the MS thesis and find the eigen values (it should be stable) with a closed loop (add a step input-a step graph should be obtained)

Step input tutorial

http://www.cambridge.org/us/features/chau/matlab/sup3.html

Matlab code for MAVs

http://www.dtic.mil/dtic/tr/fulltext/u2/a540013.pdf

Simulink tutorials

http://www.engin.umich.edu/class/ctms/pid/pid.htm

http://people.bath.ac.uk/enpms/study/docs/sem03sm_pres.pdf

Great matrix papers:

http://gust.engin.umich.edu/Papers/2010/AIAA-2010-2887.pdf

http://deepblue.lib.umich.edu/bitstream/2027.42/90638/1/AIAA-54288-189.pdf

http://www.dtic.mil/cgi-bin/GetTRDoc?AD=ADA468442

m=0.096

bx=0

by=0.12

btheta=8.0516e-4 % I is the moment of inertia for a solid cuboid

g=9.8

Kt=0.0032

I=2.6839e-4

A=[0 1 0 0 0 0; 0 -bx/m 0 0 -g 0; 0 0 0 1 0 0; 0 0 0 -by/m 0 0; 0 0 0 0 0 1; 0 0 0 0 0 -btheta/I];

B=[0 0; 0 0; 0 0; 1/m 0; 0 0; 0 -Kt/I];

C=[0 1 0 0 0 0;0 0 1 0 0 0;0 0 0 0 1 0];

K=[-390.7551 -61.4923 19.9253 2.5061 41.4851 1.0601; 441.5719 92.4694 1.2553 0.0886 -92.2723 -4.1380];

A-B*K

eig(A-B*K)

>> eig

m =

0.0960

bx =

0

by =

0.1200

btheta =

8.0516e-004

g =

9.8000

Kt =

0.0032

I =

2.6839e-004

ans =

1.0e+003 *

0 0.0010 0 0 0 0

0 0 0 0 -0.0098 0

0 0 0 0.0010 0 0

4.0704 0.6405 -0.2076 -0.0274 -0.4321 -0.0110

0 0 0 0 0 0.0010

5.2648 1.1025 0.0150 0.0011 -1.1002 -0.0523

ans =

-7.2906 + 8.4051i

-7.2906 - 8.4051i

-18.4976 + 8.8552i

-18.4976 - 8.8552i

-14.0580 + 4.4178i

-14.0580 - 4.4178i

m=0.096

bx=0

by=0.12

btheta=8.0516e-4 % I is the moment of inertia for a solid cuboid

g=9.8

Kt=0.0032

K1=2.34

I=2.6839e-4

A=[0 1 0 0 0 0; 0 -bx/m 0 0 -g 0; 0 0 0 1 0 0; 0 0 0 -by/m 0 0; 0 0 0 0 0 1; 0 0 0 0 0 -btheta/I];

B=[0 0; 0 0; 0 0; 1/m 0; 0 0; 0 -Kt/I];

C=[0 1 0 0 0 0;0 0 1 0 0 0;0 0 0 0 1 0];

K=[-390.7551 -61.4923 19.9253 2.5061 41.4851 1.0601; 441.5719 92.4694 1.2553 0.0886 -92.2723 -4.1380];

A-B*K

eig(A-B*K) %The Eigenvalue is the pole of G

t=0:.2:12;

G=tf(K1,[1 1.28]) %time constant = 1, pole at -1.28

y1=step(G,t);

plot(t,y1)

legend('pole at -1.28')

%Effect of changing the steady state gain

damp(G)

pole(G)

or

[v,d]=eig(A-B*K)

%The eigenvalues are the diagonal of the "d" matrix

%The eigenvectors are the columns of the "v" matrix

More precisely, if all eigenvalues are negative real numbers or complex numbers with negative real parts then the point is a stable attracting fixed point, and the nearby points converge to it at an exponential rate, cfLyapunov stability and exponential stability. If none of the eigenvalues is purely imaginary (or zero) then the attracting and repelling directions are related to the eigenspaces of the matrix A with eigenvalues whose real part is negative and, respectively, positive. Analogous statements are known for perturbations of more complicated orbits.

http://www.youtube.com/watch?v=i87kTpxjjsc#start=0:00;end=7:56;autoreplay=true;showoptions=false

http://www.math.siu.edu/matlab/tutorial3.pdf

http://lpsa.swarthmore.edu/MtrxVibe/EigMat/MatrixEigen.html

m1 = 2500;

m2 = 320;

k1 = 80000;

k2 = 500000;

b1 = 350;

b2 = 15020;

A=[0 1 0 0

-(b1*b2)/(m1*m2) 0 ((b1/m1)*((b1/m1)+(b1/m2)+(b2/m2)))-(k1/m1) -(b1/m1)

b2/m2 0 -((b1/m1)+(b1/m2)+(b2/m2)) 1

k2/m2 0 -((k1/m1)+(k1/m2)+(k2/m2)) 0];

B=[0 0

1/m1 (b1*b2)/(m1*m2)

0 -(b2/m2)

(1/m1)+(1/m2) -(k2/m2)];

C=[0 0 1 0];

D=[0 0];

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

Aa=[0 1 0 0 0

-(b1*b2)/(m1*m2) 0 ((b1/m1)*((b1/m1)+(b1/m2)+(b2/m2)))-(k1/m1) -(b1/m1) 0

b2/m2 0 -((b1/m1)+(b1/m2)+(b2/m2)) 1 0

k2/m2 0 -((k1/m1)+(k1/m2)+(k2/m2)) 0 0

0 0 1 0 0];

Ba=[0 0

1/m1 (b1*b2)/(m1*m2)

0 -(b2/m2)

(1/m1)+(1/m2) -(k2/m2)

0 0];

Ca=[0 0 1 0 0];

Da=[0 0];

sys=ss(Aa,Ba,Ca,Da);

K = [0 2.3e6 5e8 0 8e6]

t = 0:0.01:2;

sys_cl = ss(Aa-Ba(:,1)*K,-0.1*Ba,Ca,Da);

step(sys_cl*[0;1],t)

title('Closed-loop response to a 0.1 m step')

Introducing the reference input

Now, we will take the control system as defined above and apply a step input (we choose a small value for the step, so we remain in the region where our linearization is valid). Replace t,u and lsim in your m-file with the following,

- t = 0:0.01:2; u = 0.001*ones(size(t)); lsim(A-B*K,B,C,0,u,t)

The system does not track the step well at all; not only is the magnitude not one, but it is negative instead of positive!

Recall the schematic above, we don't compare the output to the reference; instead we measure all the states, multiply by the gain vector K, and then subtract this result from the reference. There is no reason to expect that K*x will be equal to the desired output. To eliminate this problem, we can scale the reference input to make it equal to K*x_steadystate. This scale factor is often called Nbar; it is introduced as shown in the following schematic:

We can get Nbar from Matlab by using the function rscale (place the following line of code after K = ...).

- Nbar=rscale(A,B,C,0,K)

Note that this function is not standard in Matlab. You will need to copy it to a new m-file to use it. Click here for more information on using functions in Matlab. Now, if we want to find the response of the system under state feedback with this introduction of the reference, we simply note the fact that the input is multiplied by this new factor, Nbar:

- lsim(A-B*K,B*Nbar,C,0,u,t)

and now a step can be tracked reasonably well.

Example: A State-space Controller for a Bus Suspension System

Designing the full-state feedback controller

Plotting the closed-loop response

From the main problem, the dynamic equations in state-space form are the following:

For the original problem setup and the derivation of the above equations, please refer to the Modeling page.

We want to design a feedback controller so that when the road disturbance (W) is simulated by a unit step input, the output (X1-X2) has a settling time less than 5 seconds and an overshoot less than 5%. For example, when the bus runs onto a 10 cm high step, the bus body will oscillate within a range of +/- 5 mm and will stop oscillating within 5 seconds.

The system model can be represented in MATLAB by creating a new m-file and entering the following commands (refer to main problem for the details of getting those commands). We need to define the A, B, C, D matrices by entering the following into the m-file: m1 = 2500; m2 = 320; k1 = 80000; k2 = 500000; b1 = 350; b2 = 15020; A=[0 1 0 0 -(b1*b2)/(m1*m2) 0 ((b1/m1)*((b1/m1)+(b1/m2)+(b2/m2)))-(k1/m1) -(b1/m1) b2/m2 0 -((b1/m1)+(b1/m2)+(b2/m2)) 1 k2/m2 0 -((k1/m1)+(k1/m2)+(k2/m2)) 0]; B=[0 0 1/m1 (b1*b2)/(m1*m2) 0 -(b2/m2) (1/m1)+(1/m2) -(k2/m2)]; C=[0 0 1 0]; D=[0 0]; sys=ss(A,B,C,D);

Designing the full-state feedback controller

First, let's design a full-state feedback controller for the system. Assuming for now that all the states can be measured (this assumption is probably not true but is sufficient for this problem), the schematic of the system should be:

The characteristic polynomial for this closed-loop system is the determinant of (sI-(A-B[1,0]'K)). Note that it's not sI-(A-BK) because the controller K can only control the force input u but not the road disturbance W. Recall that our B matrix is a 4 x 2 matrix, and we only need the first column of B to control u.

For this example, we have to use integral action to achieve zero steady-state error, so we add an extra state which is

. In reality the bus will eventually reach an equilibrium that yields a zero steady-state error. The new states are X1,

Y1, and Y2. Also the state-space matrices, A,B,and C, after adding extra state change to be the following:Aa=[0 1 0 0 0 -(b1*b2)/(m1*m2) 0 ((b1/m1)*((b1/m1)+(b1/m2)+(b2/m2)))-(k1/m1) -(b1/m1) 0 b2/m2 0 -((b1/m1)+(b1/m2)+(b2/m2)) 1 0 k2/m2 0 -((k1/m1)+(k1/m2)+(k2/m2)) 0 0 0 0 1 0 0]; Ba=[0 0 1/m1 (b1*b2)/(m1*m2) 0 -(b2/m2) (1/m1)+(1/m2) -(k2/m2) 0 0]; Ca=[0 0 1 0 0]; Da=[0 0]; sys=ss(Aa,Ba,Ca,Da);

Actually, there is a shortcut for MATLAB to achieve the same result.Aa = [[A,[0 0 0 0]'];[C, 0]]; Ba = [B;[0 0]]; Ca = [C,0]; Da = D; sys=ss(Aa,Ba,Ca,Da);

Add the above MATLAB code into the m-file. In this case, we treat the problem like a PID controller design. The integral control is obtained from the new state. The proportional control is obtained from a gain on Y1 or X1-X2. The direct derivative control of the output isn't possible, since derivative of Y1 or X1-X2 isn't a state. Instead we use the derivative of X1, which is available for feedback. (While X1 maybe hard to measure,

could be obtained by integrating the output of an accelerometer mounted on the bus.) It is similar to adding more damping to velocity of oscillation of the bus. Add the following MATLAB code for controller K in the m-file:K = [0 2.3e6 5e8 0 8e6]

We arrived at this value of the K, matrix by trial and error, adjusting the gain for derivative of X1,Y1 and integral of Y1, as previously mentioned.

Plotting the closed-loop response

Looking at the schematic above again, we see that after adding the K matrix into the system, the state-space equations become:

We can now obtain the closed-loop response by simply adding the following code into the m-file. Note that we need to multiply B matrix by 0.1 to simulate 0.1 m high step disturbance:t = 0:0.01:2; sys_cl = ss(Aa-Ba(:,1)*K,-0.1*Ba,Ca,Da); step(sys_cl*[0;1],t) title('Closed-loop response to a 0.1 m step')

Running the m-file in the command window, you should see the following plot:

From the plot we see that the percent overshoot and settling time requirements are satisfied. Moreover the steady-state error approaches zero as well. Therefore, we will determine that the response is satisfactory. Feel free to play around with the gain for matrix K.