%Lecture 8 on second ordersystems, bode plots, nyquist
bode(50,[1 9 30 40]);
%margin(50,[1 9 30 40]);
hold on
bode(100*50,[1 9 30 40]);
margin(50*100,[1 9 30 40]);
kp=50;num=1;den=[1 9 30 40];
[numCL,denCL]=cloop(kp*num,den);
step(numCL,denCL);
hold on
kp=100*50;num=1;den=[1 9 30 40];
[numCL2,denCL2]=cloop(kp*num,den);
step(numCL2,denCL2)
%----------------------------------------------
num=1;den=[1 1]; bode(num,den)
num=1;den=[1 10]; bode(num,den)
num=1;den= conv([1 1],[1 10]); bode(num,den);
num=2000*[1 0.5];den=conv([1 0],[1 60 500]); bode(num,den)
num=2000*[1 ];den=conv([1 0],[1 60 500]); bode(num,den)
%---------------------------------------------
%Next problem
num=10;
den=[1.25 1];
bode(num,den)

%Next
[numc,denc]=cloop(num,den)
step(numc,denc)
%--------------------------------------
%bdnawidth frequency
bode (1, [1 0.5 1 ])
%1st set of commands
w= 0.3;
	num = 1;
	den = [1 0.5 1 ];
	t=0:0.1:100; 
	u = sin(w*t);
	[y,x] = lsim(num,den,u,t);
	plot(t,y,t,u)
	axis([50,100,-2,2])

%2nd set of commands
w = 3;
	num = 1;
	den = [1 0.5 1 ];
	t=0:0.1:100; 
	u = sin(w*t);
	[y,x] = lsim(num,den,u,t);
	plot(t,y,t,u)
	axis([90, 100, -1, 1])
   
   %-------------------design---------------------------
   %Next problem
num=10;
den=[1.25 1];
bode(num,den)

%There are several several characteristics of the system that can be read directly from 
%this Bode plot. First of all, we can see that the bandwidth frequency is around 10 rad/sec. 
%Since the bandwidth frequency is roughly the same as the natural frequency (for a second
%order system of this type), the rise time is 1.8/BW=1.8/10=.18 seconds. This is a rough
%estimate, so we will say the rise time is about .2 seconds. 

%The phase margin for this system is approximately 95 degrees. This corresponds
%to a damping of PM/100=95/100=0.95. Plugging in this value into the equation 
%relating overshoot and the damping ratio (or consulting a plot of this relation),
%we find that the damping ratio corresponding to this overshoot is approximately 1%. 
%The system will be close to being overdamped. 

%The last major point of interest is steady-state error. The steady-state error can 
%be read directly off the Bode plot as well. The constant (Kp, Kv, or Ka) are located
%at the intersection of the low frequency asymptote with the w=1 line. Just extend
%the low frequency line to the w=1 line. The magnitude at this point is the constant.
%Since the Bode plot of this system is a horizontal line at low frequencies (slope = 0),
%we know this system is of type zero. Therefore, the intersection is easy to find.
%The gain is 20dB (magnitude 10). What this means is that the constant for the error
%function it 10. 
%The steady-state error is 1/(1+Kp)=1/(1+10)=0.091. If our system was type one
%instead of type zero, the constant for the steady-state error would be found
%in a manner similar to the following 

%Let's check our predictions by looking at a step response plot. This can be
%done by adding the following two lines of code into the Matlab command window. 

[numc,denc] = cloop(num,den,-1);
step(numc,denc)


%As you can see, our predictions were very good. The system has a rise time of
%about .2 seconds, is overdamped, and has a steady-state error of about 9
%. Now we need to choose a controller that will allow us to meet the design criteria. 
%We choose a PI controller because it will yield zero steady state error for a step input. 
%Also, the PI controller has a zero, which we can place. 
%This gives us additional design flexibility to help us meet our criteria. 
%Recall that a PI controller is given by: 


        K*(s+a)
Gc(s) = -------
       s

The first thing we need to find is the damping ratio corresponding to a percent overshoot
of 40%. Plugging in this value into the equation relating overshoot and damping ratio
(or consulting a plot of this relation), we find that the damping ratio corresponding 
to this overshoot is approximately 0.28. Therefore, our phase margin should be approximately 30 degrees. From our Ts*Wbw vs damping ratio plot, we find that Ts*Wbw ~ 21. We must have a bandwidth frequency greater than or equal to 12 if we want our settling time to be less than 1.75 seconds which meets the design specs. 
Now that we know our desired phase margin and bandwidth frequency, we can start our 
design. Remember that we are looking at the open-loop Bode plots. Therefore, our bandwidth
frequency will be the frequency corresponding to a gain of approximately -7 dB. 

Let's see how the integrator portion of the PI or affects our response. Change your
m-file to look like the following (this adds an integral term but no proportional term): 

num = [10];
den = [1.25, 1];
numPI = [1];
denPI = [1 0];
newnum = conv(num,numPI);
newden = conv(den,denPI);	
bode(newnum, newden, logspace(0,2))


Our phase margin and bandwidth frequency are too small. We will add gain and phase with
a zero. Let's place the zero at 1 for now and see what happens. Change your m-file to look like the following: 

num = [10];
den = [1.25, 1];
numPI = [1 3];
denPI = [1 0];
newnum = conv(num,numPI);
newden = conv(den,denPI);	
bode(newnum, newden, logspace(0,2))

It turns out that the zero at 1 with a unit gain gives us a satisfactory answer. 
Our phase margin is greater than 60 degrees (even less overshoot than expected) and
our bandwidth frequency is approximately 11 rad/s, which will give us a satisfactory response.
Although satisfactory, the response is not quite as good as we would like. Therefore,
let's try to get a higher bandwidth frequency without changing the phase margin too much.
   Let's try to increase the gain to 5 and see what happens .This will make the gain shift
      and the phase will remain the same. 

num = [10];
den = [1.25, 1];
numPI = 5*[1 1];
denPI = [1 0];
newnum = conv(num,numPI);
newden = conv(den,denPI);	
bode(newnum, newden, logspace(0,2))

a=3;T=0.0271;
cnum=[T*a 1];cden=[T 1];
num = [10];
den = [1.25, 1];
numPI = [T*a 1]/a;
denPI = [T 1 0];
newnum = conv(num,numPI);
newden = conv(den,denPI);	

That looks really good. Let's look at our step response and verify our results. 
Add the following two lines to your m-file: 

[clnum,clden] =cloop(newnum,newden,-1);
step(clnum,clden) 

%---------------------------------Bode cruise---------------------------------------

The open-loop transfer function for this problem is : 

m=1000 
b=50 
U(s)=10 
Y(s)=Velocity output 
The design criteria are: 


Rise time < 5 sec
Overshoot < 10%
Steady state error < 2%


