From: Roger Alves on
I want to include another state with the integral control (as described ln Norman Nise chapter 12), how do I do it? My file is as follows. Also is it possible to monitor the motor current through this file? How?
g=9.81; %Gravity(m/s^2)
r=0.1; %Radius of wheel(m)
Mw=0.7; %Mass of wheel(kg)
Mp=80; %Mass of body(kg)
Iw=0.0004; %Inertia of the wheel(kg*m^2)
Ip=0.04; %Inertia of the body(kg*m^2)
l=0.35; %Length to the body's centre of mass(m)
%Motor's variables
Km = 0.3035; %Motor torque constant (Nm/A)
Ke = 0.0277; %Back EMF constant (Vs/rad)
R = 1; %Nominal Terminal Resistance (Ohm)
beta = (2*Mw+(2*Iw/r^2)+Mp)
alpha = (Ip*beta + 2*Mp*l^2*(Mw + Iw/r^2))
A = [0 1 0 0;
0 (2*Km*Ke*(Mp*l*r-Ip-Mp*l^2))/(R*r^2*alpha) (Mp^2*g*l^2)/alpha 0;
0 0 0 1;
0 (2*Km*Ke*(r*beta - Mp*l))/(R*r^2*alpha) (Mp*g*l*beta)/alpha 0]
B = [ 0;
(2*Km*(Ip + Mp*l^2 - Mp*l*r))/(R*r*alpha);
0;
(2*Km*(Mp*l-r*beta)/(R*r*alpha))]
C = [1 0 0 0;
0 0 1 0]
D = [0;
0]
%TRANSFER FUNCTION FORM OF THE STATE SPACE MODEL
disp('Transfer Function of the system')
[num,den] = ss2tf(A,B,C,D)
%Obtaining the eigenvalues of the system matrix
disp('The eigenvalues of the system matrix A')
disp('A positive value will indicate an unstable system')
p = eig(A)
%%%%%%%%%%%%%%%%%
%Pole Placement control design
%%%%%%%%%%%%%%%%%
disp('The system matrix have to be full rank for pole placement control')
rank_=rank(ctrb(A,B))
P = [-9 -10 -11 -13]
K = place(A,B,P)
%%%%%%%%%%%%%%
%Simulate the system
%%%%%%%%%%%%%%
%Simulation time step
T=0:0.02:5;
%Impulse response input
U=zeros(size(T));
U(1)= 1;
%System matrices with feedback
Ac = [(A-B*K)];
Bc = [B];
Cc = [C];
Dc = [D];
%Obtaining the States and the output response
[Y,X]=lsim(Ac,Bc,Cc,Dc,U,T);
%Obtaining the torque needed to control the system
[n m] = size (X);
for i = 1:n
UU(i) = -K*X(i,:)';
end
figure
%plot the states
title('Impulse response of the plant with Pole Placement control')
plot(1), plot(T,[X(:,1) X(:,2) X(:,3) X(:,4)]), xlabel('Time [s]'),
ylabel('Position[m],Velocity[m/s],Angle[rad],Angular velocity[rad/s]')
legend('x','dx/dt','phi','dphi/dt')

Regards