Embedded Files
Download for ti-89 to convert Z(u(k-2)) = 1/(z-1)
in ti-89 u is h that represents a unit step response.
http://www.ticalc.org/archives/files/fileinfo/296/29629.html
Step response when damping ratio is one the overshoot is zero.
http://www.facstaff.bucknell.edu/mastascu/econtrolhtml/sysdyn/sysdyn2.html
http://users.tkk.fi/psipari/stadia/k05/Figs58&720.pdf
http://www.facstaff.bucknell.edu/mastascu/econtrolhtml/sysdyn/sysdyn2.html
Complicated root locus from Ogada number 3
Complicated root locus from Ogada number 2
Complicated root locus from Ogada number 1
----------------------------------------------------------------------------------------------------------------------------
SISO tool
siso tool for pole placement!
For finding that optimal gain value....
close all
clear all
gp=zpk([],[0 -2 -3],0.0707)
T=0.001
Gz=c2d(gp,T)
rlocus(Gz)
sisotool
---------------------------------------------------------------
Reduced order observers
Example 9.7.
We will again considr the desing of the system of Example 9.1.
The plant modelis given by
x(k+1)=[1 0.0952]x(k) + [ 0.00484]u(k)
[0 0.905] [0.0952]
y(k)=[1 0]x(k)
Thus we are measuring position, x1(k) and will estimate
velocity, x2(k). The closed-loop characteristic equation
was chosen to be
alphac(z) = z^2-1.776z+0.819 =0
alphae(z) = z-0.819=0
The plant state equations
Aaa=1
Aab=0.0952
Ba=0.00484
Aba=0
Abb=0.905
Bb=0.0952
Now we have the terms Ackermann's formula
G = alphae(Abb)[Aab]^-1[1]
= [0.905-0.819][0.0952]^-1[1]
= 0.903
The observer equation
q(k+1) = [0.905-(0.903)(0.0952)]q(k) + 0.903y(k+1)
+ [0-(0.903)(1)(k)+[0.0952-(0.903)(0.00484)]u(k)
q(k+1)=0.819q(k)+0.903y(k+1)-0.903y(k)+0.0908u(k)
q(k) is the estimate of velocity,x2(k).
q(k)=0.819q(k-1)+0.903y(k)-0.903y(k-1)+0.0908u(k-1)
and q(k) is the estimate at the present time. From
example 9.1 the control law is given by
u(k) = -4.52x1(k)-1.12x2(k)
which is implemented as
u(k) = -4.52y(k)-1.12q(k)
We can write the observer equation as
q(k+1) = 0.819q(k)+0.903y(k+1)-0.903Y(K)
+0.0908[-4.42y(k)-1.12q(k)]
q(k+1)=0.717q(k)+0.903y(k+1)-1.313y(k)
The control system is then implemented as follows.
A measurement y(k) is made at t=kT.
The observer state is calculated from
q(k) = 0.717q(k-1)+0.903y(k)-1.313y(k-1)
Then the control input is calculated, using
u(k) = -4.52y(k)-1.12q(k)
Example 9.8
We will now calculate the transfer function of the controller estimator
for the sytem of Example 9.7, G=0.903 and
Aaa=1
Aab=0.0952
Ba=0.00484
K1=4.52
Aba=0
Abb=0.905
Bb=0.0952
Kb=1.12
and from
Dce(z) = -U(z)/Y(z)
= K1+kb[z-Abb+GAab+(Bb-GBa)Kb]^-1
[Gz+{Aba-GAaa-K1(Bb-GBa)}]
= 4.52+1.12[z-0.905+(0.903(0.0952)+{0.052-(0.903)(0.00484)}1.12]^-1
{0.903z+[0-(0.903)(1)-4.52{0.0952-(0.903)(0.00484)}]
= 4.52+(1.12(0.903z-1.314)/(z-0.717))
= (5.53-4.71)/(z-0.717)
= 5.53(z-0.952)/(z-0.717)
Since the zero is closer to z=1 than is the pole, the control-estimator
is phase lead. The must
Example 9.9
We will consider the system of the earlier rexamples.
x(k+1) = [1 0.0952]x(k)+[0.00484]u(k)
[0 .0905] [0.0952]
y(k) = [1 0]x(k)
with optimal gains k =[4.52 1.12]
we will use the same observer characteristic equation as
used for the prediction observer of Example 9.3
alphae(z) = z^2-1.638z+0.671
From example 9.3,
alphae(A) = [0.033 0.0254]
[0 0.00763]
also
CA=[1 0.0952]
and thus
CA^2 = (CA)A = [1 0.0953][1 0.0952]
[0 0.905]
= [1 0.1814]
Then
[CA]^-1 = [1 0.0952]^-1= [ 2.104 -1.104]
[CA^2] [1 0.1814] [-11.60 11.60]
From (9-72)
G=alphae(A)[CA]^-1[0]
[CA^2][1]
= [0.033 0.0254][2.104 -1.104][0]
[0 0.00763 ][-11.60 11.60][1]
= [0.033 0.0254][-1.104]
[0 0.00763 ][11.60]
= [0.258]
[0.0885]
The difference equation of the observer in (9-68) will now be evaluated
GCA=[0.258][1 0][1 0.0952]
[0.0885][0 0.905 ]
= [0.258][1 0.0952]
[0.0885]
= [0.258 0.0246]
[0.0885 0.00843]
A-GCA = [0.742 0.0706 ]
[-0.0885 0.897]
Also, evaluation of (B-GCB) yields
B-GCB = [0.00359]
[ 0.0948]
we see then from 9-68 that the observer is implemented as
q(k+1) = [0.742 0.0706]q(k) + [0.00359]u(k) + [0.258]y(k+1)
[-0.0885 0.897] [0.0948 ] [0.0885]
Example 9.10.
This example is a continuation of the last exampl. The transfer function of the controller
estimator will be calculated from (9-73) from Example 9.9.
A-GCA = [0.742 0.0706]
[-0.0885 0.897]
B-GCB= [0.00359]
[0.0948 ]
hence
(B-GCB)K = [0.00359][4.52 1.12]
[0.0948]
= [0.0219 0.00542]
[0.4303 0.1066]
and
A-GCA-BK+GCBK =[0.726 0.0666]
[-0.517 0.791]
[zI-A+GCA+BK-GCBK]^-1=[z-0.726 -0.0666]^-1
[ 0.157 z-0.791]
=(1/delta)[z-0.791 0.0666]^-1
[-0.157 z-0.726]
delta=[zI-A+GCA+BK-GCBK]
=z^2-1.52z+0.609
Then from (9-73)
Dce(z) = zK[zI-A+GCA+BK-GCBK]^-1*G
=(z/delta)[4.52 1.12][z-0.791 0.0666][0.258]
[-0.517 z-0.726][0.0885]
=(z/delta)[4.52 1.12][0.258z - 0.198]
[0.0885z-0.1976]
=(1.27z^2-1.12z)/(z^2-1.52z+0.609)
A problem can occur in the implementation of obsever-based controllers. Thus
far we have not considered the relative stability of these control system.
A system that has adequate stability margins is said to be robust; obsever-
based control system may not be roboust in terms of the gain and phase
margins that appear at the input to the plant. To determine these stability margins,
we see from Figure 9-8 that the frequency response for the open-loop
function Dce(z)G(z) must be calculated. This was done for the prediction-observer system.
Kalman filters are optimal currenet observers.
The state model of the closed-loop system employing the current
observer can be derived as
[x(k+1)] = [A -BK][x(k)]
[q(k+1)] = [GCA A-GCA-BK][q(k)]
The derivation of this equation is straightforward and is given
as an exercise in Problem 9-30.
--------------------------------------------------------------------------
hw7_digitalcontrols
Problem1.
Assume that a closed-loop characterisitcs equation is given by
(z-p1)(z-p2)...(z-pn)=0. Show that the requirement for
system stability is that the magnitudes of all poles
of the closed-loop transfer function are less than unity,
that is |pi|<1, i=1,2...n.
c(z)=k1z/(z-p1) + ...knZ/(z-pn)
z^-1(c(z)) = z^-1(k1z/(z-p1)) +...z^-1(knz/(z-pn))
c(k) = k1p1^k + ...knPn^k
limk>infinity k1P1^k is bounded for |P1|<1
Problem 2.
Consider a sampled data system with T=0.5s and the
characteristics equation given by:
(z-0.9)(z-0.8)(z^2-1.9z+1)=0
a)Find the poles of the system
b) A discrete linear-time invariant system is stable,
unstable, or margianlly stable. Identify the type
of stability for this system.
c) The natural response of this system contains an
undamped sinusoidal response term of the form
Acos(wkT+theta). Find the frequency w of this term.
a) z=0.9, z=0.8, TI-89> zeros(x^2-1.9x+1,x)=0
z=0.95+-0.312j
r=sqrt(0.95^2+0.3122^2) = 1
theta = tan^-1(0.3122/0.95)
= +-18.19degree
18.19deg(pi/180) = 0.3175rad/sec
b) z=0.8,z=0.9
z=r<theta = 1<+-18.19degree
marginally stable
c) z = r<theta
= r<wT
wT=0.3175rad/sec
T=0.5sec
w = 0.3175/0.5
= 0.635 rad/sec
w = 0.635rad/sec
Problem 3.
Consider athe system of figure 1 and let the digital controller
be a vairable gain K, such that D(z) = K.
a) Write the closed-loop system characteristic equation as a
function of the sample period T
R
__0__/ __D(z)>>(1-e^-ts/s)---(1/s+1)---c
| |
-------------------------------------
a) G(z)/1+kG(z) = 0
1+kG(z) = 0 is the closed loop characteristic equation
G(z) = (z-1)/z*z(Gp(s)/s)
= (z-1)/z Z((1/s)(1/(s+1)))
PFE, TI-89 expand(1/x(1/x+1)) = 1/x -1/(x+1)
G(z) = (z-1)/z* Z[1/s - 1/(s+1)]
= (z-1)/z*[z/(z-1) - z/(z-e^-T]
= 1-(z-1)/(z-e^-T)
= e^-T+1/(z-e^-T)
1+kG(z) = 0
1+k(-e^-T+1/(z-e^-T)) =0
z-e^-T-ke^-T+k=0
z-e^-T+k(-e^-T+1)=0
b) Determine the ranges of K>0 for stability for the sample periods
T=1s, T=0.1s, T=0.01s.
characteristic equation:
z-e^-T+k(-e^-T+1)=0
for T=1sec
>> z-0.3678+k(0.632)=0
z=0.33678-0.632k >-1
k<2.164
for T=0.1sec
>> z-0.9048+k(0.0951)=0
z=0.9048-0.0951k > -1
k<20.02
for T= 0.01 sec
>> z-0.99005+k(0.00995)=0
z=0.9905-0.00995k >-1
k<200.05
c) Consider the system with all sampling removed and with
Gp(s) = k/(s+1). Find the range of K>0 for which the analog
system is stable.
Gp(s)/(1+kGp(s))=0
closed loop characterisitc equation
1+kGp(s) = 0
1+k/(s+1) = 0
s+1+k=0
s=-k-1
if k>-1 system is stable
for k>0, the rang of k is from
(0,infinity)
d) Comparing the ranges of K from b and c, give the effects on
stability of reducing the sample period T.
By reducing the sampling period T will increase
the gain, k and make the system more stable!
T dec, k inc, system becomes more stable.
Problem 4.
Consider the robot arm joint controls system
of Fig 2. T =0.1s and D(z) =1.
Let Z[(1-e^-Ts)/s(4/s(s+2)] = 0.01873z+0.01752/(z-1)(z-0.8187)
thetac
Controller
__0__/ __D(z)>>ZOH>>K>>(200/(s(0.5s+1))---(1/100)---thetaa
| |
--------------volts----0.07 sensor-----------------------
a) Write the closed-loop system characteristic equation
G(z)/(1+kG(z))=0
1+kG(z)=0
1+(0.07)K((0.01873z+0.01752/(z-1)(z-0.8187))=0
(z-1)(z-0.8187)+(0.07)k(0.01873z+0.01752)=0
z^2-1.8187z+0.8187+0.00131kz+0.001226k=0
z^2-(1.8187z-0.001311k)z+0.08187+0.001226k=0
This is the closed-loop characteristic equation.
b) use the Routh-Hurwitz criterion to determine the range
of K for stability
z=1+(T/2)w/(1-(T/2)w), T=0.1sec
z=(1+0.05w)/(1-0.05w)
= 20+w/(20-w)
z^2-(1.8187-0.001311k)z+0.8187+0.001226k=0
Q(w) = (20+w)^2-(1.8187-0.001311k)(400-w^2)
w^2 | 3.6374-0.000085k >0 k<42,793
w^1 | 7.252-0.04904k>0 K<147.9
w^0 | 1.0148k >0 k>0
for stability 0<k<147.9
c) Check the restuls of part b) using the Jury test.
1.
Q(1)>0
characteristic equation
z^2-(1.8187-0.001311k)z+0.8187+0.001226k=0
1^2-1.8187+0.001311k+0.8187+0.001226k>0
k>0
2. (-1)^2Q(-1)
1^2-1.8187+0.001311k+0.8187+0.001266k>0
k<42.793
3. |ao|<|a2|
-0.8187+0.001226k<1
k<147.9
to be stable
0<k<147.9
d) Determine the location of all roots of the characteristic
equation in both the w-plane and the z-plane for the value
of K>0 for which the system is marginally stable.
found k=147.9
a) z-plane plug k into characteristic equation
z^2-(1.8187-0.001311k)z+0.8187+0.001226k=0
z^2-1.6248z+1=0
TI-89 czeros(x^2-1.6248x+1,x)
>> (0.8124 +- 0.5813j)
r = sqrt(0.8124^2+0.5913^2) = 1
theta = tan^-1(0.5813/0.8124)
= +-35.67degree
z = r<theta
= 1<+-35.67degree
w-plane, k=147.9
Q(w) = (3.6374-0.000085k)w^2+(7.252-0.04904k)w+1.0148k=0
= 3.6148w^2+150.09=0
w=+-6.43j
e) Determine both the z-plane frequency and the w-plane
frequency at which the system will oscillate when
margianlly stable, using the results of part d.
z-plane
wz=r<theta
= 1<+-35.67degree
35.67degree(pi/180degree) = 0.6225 rad
<theta = wT
0.6225 = w(0.1)
w-plane
wz=6.225rad/sec
Ww=+-6.43j
Ww= 6.43rad/sec
Problem 5.
For a system described by G(z) = 0.6321/(z-0.3679)
a) Plot (by hand) the z-plane root locus
1pole-0zeros=1 go to infinity
Imaginary
|
|
<<<<<--x---Real
| 0.3679
|
|
b) Plot (by hand) the w-plane root locus
z=(1+(T/2)w)/(1-(T/2)w)
= 2+w/(2-w)
G(w) = 0.6321(2-w)/(2+w)-0.3679(2-w)
= -0.4621(w-2)/(w+0.9242)
Imaginary
| w-plane root locus
|
<<<x<<-|---0<<<<-Real
-0.924 | 2
|
|
c) Determine the range of K for stability using
the results of part a.
k for z-plane , z=-1
k = abs(1/g(z))
= abs(z-0.3679)/0.6321
= abs(-1-0.3679/0.6321)
= abs(-2.164)
= 2.164
d) Determine the range of K for stability using
the results of part b.
k for w-plane , z=infinity
k = abs(1/g(w))
= abs(w+0.9242)/(-0.4621(w-2))
= abs(infinity+0.9242/(-0.4621(w-2))
= abs(1/-0.462)
= 2.164
Cancel poles and add them zeros
z transform of 1/s^3
http://control.dii.unisi.it/sdc/altro/TabellaTrasformataZ.pdf
Determining Stability:
Signal Flow diagram
-----------------------------------------------------------------------------------------------------------------------
Matlab Notes for digital controls class:
%
% 2/11/14, lecture notes
% digital controls
%
clear all
close all
for k=0:1:100 %a lot of samples
c(k+1)=0.5*(1-(1-0.2642).^k)
%this is the discrete system
% difference equation
figure(1)
plot(c)
grid minor
figure(2)
stem(c)
end
grid minor
xlabel('number of samples')
title('Response for D(Z)=1')
ylabel('Amplitude')
for n=0:1:10
T(n+1)=0.5*(1-(1-0.2642).^n)
figure(3)
stem(T)
end
grid minor
xlabel('number of samples')
title('Response for D(Z)=1')
ylabel('Amplitude')
% axis equal
% dcgain(c) % or final value
% damping(c)
-----------------------------------------------------------------------------------------------------------------------
Konstantinos Tsakalis
http://tsakalis.faculty.asu.edu/
Matlab script for generating root locus of different discrete systems:
Matlab script used to design and simulate a discrete controller
clear
% Matlab script used to design and simulate a discrete controller
% for the analog plant Gp_s defined in the Ch5 lecture notes for in-class
% example #5 with the desired time domain specifications such that
% the desired closed-loop poles are located at:
s_cl= [-5+j*5, -5-j*5]
Gp_s=tf(1,conv([1 3],[1 1])) % Define the analog plant (note its a Type 0 system w/o compensation)
figure;step(Gp_s) % Compute the unit step response of the open-loop plant (note the response is slow w/finite ess~=0)
% Use the sisotool(Gp_s) command to design a cascade PID controller w/required gain
% We can begin by making an educated guess at the form of the controller:
K=1;Gc_s=tf(K*conv([1 1],[1 7.14]),[1 0]) % Define a cascade PID controller w/unity gain
% This allows us to start the design process using the command sisotool(Gp_s,Gc_s), if desired
% Note, the angle (for calculation of the poles/zeros) and magnitude condition for calculation of the required gain)
% could have been used to do the same by design by hand, or alternately coefficient matching of the desired
% and actual characteristic polynominals
K=7.0;Gc_s=tf(K*conv([1 1],[1 7.14]),[1 0]) % Define a cascade PID controller w/required gain from RL plot
L_s=Gc_s*Gp_s % Compute the analog loop gain
zpk(L_s) % View loop gain poles and zeros to ensure no unstable cancellation and avoid internal instability
figure;plot(real(s_cl),imag(s_cl),'rs');axis([-8 1 -4.5 4.5]) % Plot the desired closed-loop poles in the s-plane
hold;rlocus(L_s) % Create RL plot of the loop gain (why does the plot intersect the desired closed-loop poles at unity gain?)
Hcl_s=feedback(L_s,1,-1) % Compute the closed-loop transfer function for the analog system design
kp=evalfr(L_s,0) % Compute the position error constant kp=L(0)
ess=1/(1+kp) % Compute ess to a step input (does this agree with the plot?)
% If the step response of the closed-loop analog system is acceptable,
% we must then convert the analog design over to the discrete domain, as shown:
T=0.02;GZAS_z=c2d(Gp_s,T,'zoh') % Discretize the analog plant w/ZOH and sampler
Gc_z=c2d(Gc_s,T,'tustin') % Discretize the analog controller using the bilinear transform
% Note, how can this be written as a difference eqn. involving e(k) and u(k)?
L_z=Gc_z*GZAS_z % Compute the discrete loop gain
zpk(L_z) % View loop gain poles and zeros to ensure no unstable cancellation and avoid internal instability
z_cl=exp(s_cl*T) % Compute the desired closed-loop poles in the z-plane
figure;plot(real(z_cl),imag(z_cl),'rs');axis([-2 2 -2 2]) % Plot the desired closed-loop poles in the s-plane
hold;rlocus(L_z) % Create RL plot of the loop gain (why does the plot intersect the desired closed-loop poles at ~unity gain?)
% Note: You may need to zoom in on the region of the RL plot around the z_cl to answer the question posed above
Hcl_z=feedback(L_z,1,-1) % Compute the closed-loop transfer function for the discrete system design
figure;step(Hcl_s);hold;step(Hcl_z); % Compare the closed-loop step responses of the analog and discrete systems
% Repeat the process until an acceptable design is achieved, if necessary
Coomet Annimation:
%G=zpk([1 10],[-1 2 -5 0.5],1)
G=zpk([1],[0.5 -1+j -1-j -3],3)
theta=[-pi:pi/2000:pi]';
%s1=(2+0.25*cos(6*theta)).*exp(-j*theta);
s1=1+0.25*cos(-theta)+j*sin(-theta);
s2=-1+j+(0.8+0.1*sin(3*theta)).*exp(-j*theta);
s3=-1+3*cos(-theta)+2*j*sin(-theta);
%s3=j*10*theta;
for ii=1:size(theta,1)
F1(ii)=evalfr(G,s1(ii));
F2(ii)=evalfr(G,s2(ii));
F3(ii)=evalfr(G,s3(ii));
s4(ii)=5*(max(0,cos(-theta(ii)))+j*sin(-theta(ii)));
F4(ii)=evalfr(G,s4(ii));
end
figure(1)
hold off,
plot(s1),grid,
hold on,plot(s2,'r'),plot(s3,'k'),plot(s4,'g'),
pzmap(G)
axis equal
pause
figure(2),
hold off,plot(s1),axis([-6 7 -5 5])
grid
hold on,comet(real(s1),imag(s1))
figure(3)
hold off,plot(F1),grid,
hold on,comet(real(F1),imag(F1))
pause
figure(2)
hold on, plot(s2,'r'),comet(real(s2),imag(s2))
figure(3)
hold on, plot(F2,'r'),comet(real(F2),imag(F2))
pause
figure(2)
hold on, plot(s3,'k'),comet(real(s3),imag(s3))
figure(3)
hold on, plot(F3,'k'),comet(real(F3),imag(F3))
pause
figure(2)
hold on, plot(s4,'g'),comet(real(s4),imag(s4))
figure(3)
hold on, plot(F4,'g'),comet(real(F4),imag(F4))
break
pause
G=zpk([-1],[0 3],5)
for ii=1:size(theta,1)
s4(ii)=0.01+10*(max([0,cos(-theta(ii)),0.01*cos(10*(pi-theta(ii)))])+j*sin(-theta(ii)));
F4(ii)=evalfr(G,s4(ii));
end
figure, plot(s4),grid, hold on, pzmap (G),comet(real(s4),imag(s4))
figure, plot(F4),grid, hold on,comet(real(F4),imag(F4))
------------------------------------------------------------------------------
Root Locus:
clear
close all
tf=8;
t=0:tf/250:tf;
%---phase lead compensation
locus=figure
plant=zpk([],[0,-1],3); %---plant
lead=zpk([-3],[-7],10/3) %---lag compensator
gh=series(lead,plant);
rlocus(gh)
a=axis; a(3)=-3; a(4)=3; axis(a); a=axis;
axis equal
hold on
plot([-1,-1],[-20,20],'r--')
plot([0,-20*0.4],[0,20*sin(acos(0.4))],'r--')
plot([0,-20*0.4],[0,-20*sin(acos(0.4))],'r--')
pause
r=rlocus(gh,1)
title('lead compensator')
plot(r,'gs', 'MarkerEdgeColor','k',...
'MarkerFaceColor','g',...
'MarkerSize',9)
pause
stepr=figure;
sys=feedback(gh,1);
step(sys,t);
title('lead compensator')
hold on
plot([0,8],[1.02,1.02],'r--');
plot([0,8],[.98,.98],'r--');
plot([4,4],[0,10],'r--');
plot([0,8],[1.25,1.25],'r--');
%----add phase lag compensation
pause
figure
lag=zpk([-.15],[-.015],1) %---lag compensator
gh=series(lag,gh); %---in series with gh
rlocus(gh); axis equal
title('lead and lag compensators')
axis(a)
hold on
plot([-1,-1],[-20,20],'r--')
plot([0,-20*0.4],[0,20*sin(acos(0.4))],'r--')
plot([0,-20*0.4],[0,-20*sin(acos(0.4))],'r--')
pause
r=rlocus(gh,1)
plot(r,'gs', 'MarkerEdgeColor','k',...
'MarkerFaceColor','g',...
'MarkerSize',9)
pause
figure
sys2=feedback(gh,1);
sysI=series(sys,zpk([],[0],1)); %---integrate output
sys2I=series(sys2,zpk([],[0],1)); %---integrate output
Tf=20; %---final time
step(sysI,sys2I,Tf); %---compare ramp responses
title('lead-lag compensator')
a=axis; a(2)=Tf; a(4)=Tf; axis(a);
grid
hold on
plot([0,Tf],[0,Tf],'r--')
legend('lead','lead-lag')
pause
figure %---did we damage the step response?
step(sys,sys2,t)
title('lead and lag compensators')
hold on
plot([0,8],[1.02,1.02],'r--');
plot([0,8],[.98,.98],'r--');
plot([4,4],[0,10],'r--');
plot([0,8],[1.25,1.25],'r--');
legend('lead','lead-lag')
orient landscape
print(locus, '-dpdf','rootlocus.pdf')
print(stepr, '-dpdf','step.pdf')
Aliasing Problem:
%%
fo = 4; %frequency of the sine wave
Fs = 50; %sampling rate
Ts = 1/Fs; %sampling time interval
t = 0:Ts:1-Ts;
n = length(t); %number of samples
y = 2*sin(2*pi*fo*t);
plot(t,y)
hold on
stem(t,y)
xlabel('time (seconds)')
ylabel('y(t)')
title('Sample Sine Wave')
legend('Continuous sin wave','Discrete samples')
%%
%=================================================================
%
% aliasing_example.m: the MATLAB code for illustrating aliasing
% For further explanations, refer to the course notes, "aliasing" section
% input: c -> amplitude of the sinusoid
% w_0 -> angular frequency of the sinosoid
% T_s -> sampling period
% k -> shifting amount of the new sinosoid
% output: 3 plots; one for the original sinusoid,
% one for the new sinusoid
% one for both of them
%
% This code is generated as illustration of the aliasing problem in
% the course notes.
%
%==================================================================
c=2
w_0=1
T_s=1
k=1
t = 0:0.01:10; % Get a time interval between [0,10] with 0.01 interval length
f_t = c .* exp( i * w_0 .* t ); % Calculate f(t)
nTs = 0:T_s:10; % nTs stands for n*T_s in the notes, I merged them for readability
f_n = c .* exp( i * w_0 .* nTs); % Calculate f(n)
f_0 = w_0/(2*pi); % The frequency of the original sinusoid
f_s = 1/T_s; % Sampling frequency
% g is the new sinusoid. Calculate g(t) and g(n)
g_t = c .* exp( i * 2*pi*(f_0+k*f_s) .* t );
g_n = c .* exp( i * 2*pi*(f_0+k*f_s) .* nTs);
% NOTE: In the following code, LaTeX format is used in some strings.
% First Plot: Plot the original sinusoid with a blue line and the sample
% points with blue squares.
subplot(3,1,1),
plot(t,real(f_t),'-b');
hold on, plot(nTs,real(f_n),'LineStyle','none',...
'Marker','o',...
'MarkerEdgeColor','k',...
'MarkerFaceColor','b',...
'MarkerSize',6);
title(['c = ' num2str(c) ', w_0 = \pi/' num2str(pi/w_0) ', T_s = ' num2str(T_s),...
', k = ' num2str(k), '\newline \newline Original sinusoid and its samples']);
legend('f(t) = ce^{iw_0t}','f(n) = ce^{iw_0nT_s}');
% Second Plot: Plot the new sinusoid with a red line and the sample
% points with red squares.
subplot(3,1,2),
plot(t,real(g_t),'-r');
hold on, plot(nTs,real(g_n),'LineStyle','none',...
'Marker','o',...
'MarkerEdgeColor','k',...
'MarkerFaceColor','r',...
'MarkerSize',6);
title('Another sinusoid and its samples');
legend('g(t) = ce^{i2\pi(f_0+kf_s)t}','g(n) = ce^{i2\pi(f_0+kf_s)nT_s}');
% Last Plot: Plot both sinusoids the sample points with green squares.
subplot(3,1,3),
plot(t,real(f_t),'-b');
hold on, plot(t,real(g_t),'-r');
plot(nTs,real(g_n),'LineStyle','none',...
'Marker','o',...
'MarkerEdgeColor','k',...
'MarkerFaceColor','g',...
'MarkerSize',6);
title('The two sinusoids and samples on the same plot');
legend('f(t)','g(t)','f(n) or g(n)');
ZPK to represent a discrete model for zeros,poles and gains:
clear
clc
% 2nd order - \zeta = \sqrt{2}
sys1 = zpk([],[-10+10*j -10-10*j],1);
sys2 = zpk([-10],[-10+10*j -10-10*j],1);
sys3 = zpk([],[-1+10*j -1-10*j],1);
% Butterworth filter
om_c = 10;
sys4 = zpk([],[-om_c -om_c*exp(j*pi/3) -om_c*exp(-j*pi/3)],om_c^3);
sys5 = zpk([10],[-10+10*j -10-10*j],1);
sys6 = zpk([0],[-10 -10],1);
figure(1)
subplot(321), bode(sys1,{1,100})
title('Bode plot 1'); axis([1 10^2 -80 10]);
subplot(322), bode(sys2,{1,100});
title('Bode plot 2'); axis([1 10^2 -80 10]);
subplot(323), bode(sys3,{1,100});
title('Bode plot 3'); axis([1 10^2 -80 10]);
subplot(324), bode(sys4,{1,100});
title('Bode plot 4'); axis([1 10^2 -80 10]);
subplot(325), bode(sys5,{1,100});
title('Bode plot 5'); axis([1 10^2 -80 10]);
subplot(326), bode(sys6,{1,100});
title('Bode plot 6'); axis([1 10^2 -80 10]);
figure(2)
subplot(321),pzmap(sys4); grid;
axis([-11 11 -11 11]);
title('Pole-zero map 1');grid;
subplot(322),pzmap(sys2);grid;
axis([-11 11 -11 11]);
title('Pole-zero map 2'); grid;
subplot(323),pzmap(sys1);grid;
axis([-11 11 -11 11]);
title('Pole-zero map 3'); grid;
subplot(324),pzmap(sys6);grid;
axis([-11 11 -11 11]);
title('Pole-zero map 4');grid;
subplot(325),pzmap(sys5);grid;
axis([-11 11 -11 11]);
title('Pole-zero map 5');grid;
subplot(326),pzmap(sys2);grid;
axis([-11 11 -11 11]);
title('Pole-zero map 6'); grid;
figure(3)
subplot(321),step(sys3); grid;
title('Step response 1');
subplot(322),step(sys4); grid;
title('Step response 2');
subplot(323),step(sys1); grid;
title('Step response 3');
subplot(324),step(sys6); grid;
title('Step response 4');
subplot(325),step(sys2); grid;
title('Step response 5');
subplot(326),step(sys5); grid;
title('Step response 6');
Binary numbers:
http://en.wikipedia.org/wiki/Binary_number
http://l3d.cs.colorado.edu/courses/CSCI1200-96/binary.html
0000 = 0
0001 = 1
0010 = 2
0011 = 3
0100 = 4
0101 = 5
0110 = 6
0111 = 7
1000 = 8
1001 = 9
1010 = 10
1011 = 11
1100 = 12
1101 = 13
1110 = 14
1111 = 15
Resources Bought:
k you for shopping with us. We thought you'd like to know that we shipped your items, and that this completes your order. Your order is on its way, and can no longer be changed. If you need to return an item from this shipment or manage other orders, please visit Your Orders on Amazon.com.
Your estimated delivery date is:
Friday, January 17, 2014
Your order was sent to:
Michael Thompson
2001 N 39TH PL
PHOENIX, ARIZONA 85008-3020
United States
Depending on the ship speed you chose, it may take 24 hours for tracking information to be available in your account.
Shipment Details
Item Subtotal:
Shipping & Handling:
Promotion Applied:
Total Before Tax:
Sales Tax Collected:
Shipment Total:
Paid by Visa:
$144.37
$12.00
-$0.00
$156.37
$1.80
$158.17
$158.17
Hello Michael,
Hello Michael,
Thank you for shopping with us. We'd like to let you know that Amazon has received your order, and is preparing it for shipment. Your estimated delivery date is below. If you would like to view the status of your order or make any changes to it, please visit Your Orders on Amazon.com.
Your estimated delivery date is:
Tuesday, January 21, 2014
Your shipping speed:
Standard Shipping
Your order will be sent to:
Michael Thompson
2001 N 39TH PL
PHOENIX, ARIZONA 85008-3020
United States
Line Follow, PID Controller in Arduino
Convert Continuous transfer function to Discrete transfer function
http://www.mathworks.com/help/control/ref/c2d.html
cretize the continuous-time transfer function:
with input delay Td = 0.35 second. To discretize this system using the triangle (first-order hold) approximation with sample time Ts = 0.1 second, type
H = tf([1 -1], [1 4 5], 'inputdelay', 0.35); Hd = c2d(H, 0.1, 'foh'); % discretize with FOH method and % 0.1 second sample time Transfer function: 0.0115 z^3 + 0.0456 z^2 - 0.0562 z - 0.009104 --------------------------------------------- z^6 - 1.629 z^5 + 0.6703 z^4 Sampling time: 0.1
The next command compares the continuous and discretized step responses.
step(H,'-',Hd,'--')
Discretize the delayed transfer function
using zero-order hold on the input, and a 10-Hz sampling rate.
h = tf(10,[1 3 10],'iodelay',0.25); % create transfer function hd = c2d(h, 0.1) % zoh is the default method
These commands produce the discrete-time transfer function
Transfer function: 0.01187 z^2 + 0.06408 z + 0.009721 z^(-3) * ---------------------------------- z^2 - 1.655 z + 0.7408 Sampling time: 0.1
In this example, the discretized model hd has a delay of three sampling periods. The discretization algorithm absorbs the residual half-period delay into the coefficients of hd.
Compare the step responses of the continuous and discretized models using
step(h,'--',hd,'-')
Discretize a state-space model with time delay, using a Thiran filter to model fractional delays:
sys = ss(tf([1, 2], [1, 4, 2])); % create a state-space model sys.InputDelay = 2.7 % add input delay
This command creates a continuous-time state-space model with two states, as the output shows:
a = x1 x2 x1 -4 -2 x2 1 0 b = u1 x1 2 x2 0 c = x1 x2 y1 0.5 1 d = u1 y1 0 Input delays (listed by channel): 2.7 Continuous-time model.
Use c2dOptions to create a set of discretization options, and discretize the model. This example uses the Tustin discretization method.
opt = c2dOptions('Method', 'tustin', 'FractDelayApproxOrder', 3); sysd1 = c2d(sys, 1, opt) % 1s sampling time
These commands yield the result
a = x1 x2 x3 x4 x5 x1 -0.4286 -0.5714 -0.00265 0.06954 2.286 x2 0.2857 0.7143 -0.001325 0.03477 1.143 x3 0 0 -0.2432 0.1449 -0.1153 x4 0 0 0.25 0 0 x5 0 0 0 0.125 0 b = u1 x1 0.002058 x2 0.001029 x3 8 x4 0 x5 0 c = x1 x2 x3 x4 x5 y1 0.2857 0.7143 -0.001325 0.03477 1.143 d = u1 y1 0.001029 Sampling time: 1 Discrete-time model.
The discretized model now contains three additional states x3, x4, and x5 corresponding to a third-order Thiran filter. Since the time delay divided by the sampling time is 2.7, the third-order Thiran filter (FractDelayApproxOrder = 3) can approximate the entire time delay.
Discretize an identified, continuous-time transfer function and compare its performance against a directly estimated discrete-time model
Estimate a continuous-time transfer function and discretize it.
load iddata1 sys1c = tfest(z1, 2); sys1d = c2d(sys1c, 0.1, 'zoh');
Estimate a second order discrete-time transfer function.
sys2d = tfest(z1, 2, 'Ts', 0.1);
Compare the two models.
compare(z1, sys1d, sys2d)
The two systems are virtually identical.
Discretize an identified state-space model to build a one-step ahead predictor of its response.
load iddata2 sysc = ssest(z2, 4); sysd = c2d(sysc, 0.1, 'zoh'); [A,B,C,D,K] = idssdata(sysd); Predictor = ss(A-K*C, [K B-K*D], C, [0 D], 0.1);
zplane
zplane
Zero-pole plot
http://www.mathworks.com/help/simulink/slref/discretezeropole.html
http://www.mathworks.com/help/signal/ref/zplane.html
I am currently taking digital control: design and implementation this Spring 2014 that is being
taught by professor Panagiotis Artemiadis.
Asst Professor
Sch Engr Matter Trnsprt Energy
Faculty
Mail Code: 6106
ERC 355 (Map)
(480)965-418
The following are my notes for this course:
Page updated
Google Sites
Report abuse