Menu
▾
▴
python-control-discuss — Discussion list for Python Control Systems Library
You can subscribe to this list here.
| 2010 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
(19) |
Jul
|
Aug
|
Sep
|
Oct
|
Nov
(11) |
Dec
(5) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2011 |
Jan
|
Feb
(113) |
Mar
(12) |
Apr
(27) |
May
(2) |
Jun
(16) |
Jul
(6) |
Aug
(6) |
Sep
|
Oct
(3) |
Nov
(9) |
Dec
(2) |
| 2012 |
Jan
(5) |
Feb
(11) |
Mar
|
Apr
(3) |
May
|
Jun
|
Jul
|
Aug
(3) |
Sep
(7) |
Oct
(18) |
Nov
(18) |
Dec
|
| 2013 |
Jan
(4) |
Feb
(1) |
Mar
(3) |
Apr
(1) |
May
|
Jun
(33) |
Jul
(2) |
Aug
(5) |
Sep
|
Oct
|
Nov
|
Dec
(1) |
| 2014 |
Jan
(1) |
Feb
|
Mar
(8) |
Apr
|
May
(3) |
Jun
(3) |
Jul
(9) |
Aug
(5) |
Sep
(6) |
Oct
|
Nov
|
Dec
|
| 2015 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
(4) |
| 2017 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
(5) |
Nov
|
Dec
|
| 2018 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
(1) |
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
|
| 2019 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
(6) |
Sep
|
Oct
|
Nov
(2) |
Dec
|
| 2020 |
Jan
(1) |
Feb
(1) |
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
|
| 2021 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
(4) |
Oct
|
Nov
|
Dec
|
|
From: Richard M. <mu...@cd...> - 2010-11-06 15:25:47
|
Thanks, Ryan. I'll pull the latest version of your package into the sourceforge repository (so that developers can access it quickly) and update the python-control step response to use your code. -richard On 6 Nov 2010, at 8:22 , Ryan Krauss wrote: > FYI, I saw some mention of lsim2 in a recent email. My code includes > a step response function that first tries lsim then uses lsim2 if > there is an exception. The exceptions usually are caused by pure > integrators. I think that defaulting to lsim2 will make the > simulation slower. > > I have fixed a few bugs in my code. There may even be a bug fix in > the digital conversion code that hasn't made it into my main branch. > The most recent versions of my code are on github and on my course > website for this semester: > > https://github.com/ryanGT/controls > http://www.cs.siue.edu/~rkrauss/mechatronics/2010/python_controls/python.html > > -- > Ryan Krauss, Ph.D. > Assistant Professor > Mechanical Engineering > Southern Illinois University Edwardsville > > > > On Sat, Nov 6, 2010 at 8:45 AM, Richard Murray <mu...@cd...> wrote: >> A couple of updates on status of the the python-control package, related to the commit earlier today: >> >> * Ryan Krauss has sent me his control.py package that he uses for his class at Southern Illinois University Edwardsville. He has implemented functions for continuous to discrete time conversion and root locus plus. Ryan's code is now in trunk/external/control.py and Ryan has graciously given us permission to incorporate functions into the python-control package. >> >> * Roberto Bucher has sent me his yottalab.py package that has a number of functions for state space computations in continuous and discrete time. This code is in trunk/external/yottalab.py. Roberto has also give us permission to incoporate functions into the python-control library. >> >> * I've added a pade() command developed by Sawyer Fuller at Caltech, using the description in Golub and Van Loan. This code is in trunk/src/delay.py. >> >> * A group of students at Princeton (Lauren Padilla, Kevin Chen and Steve Brunton) are working on a class project this term and will be helping to implement some additional functions in the library. I've encouraged them to use this list to send out thoughts and ideas about their plans. >> >> * In integrating some of the functions from Ryan and Roberto, I've made various changes to things and fixed a couple of bugs. See trunk/ChangeLog for details. >> >> That's it for now. Comments and feedback is always welcome. >> >> -richard >> >> >> >> ------------------------------------------------------------------------------ >> The Next 800 Companies to Lead America's Growth: New Video Whitepaper >> David G. Thomson, author of the best-selling book "Blueprint to a >> Billion" shares his insights and actions to help propel your >> business during the next growth cycle. Listen Now! >> http://p.sf.net/sfu/SAP-dev2dev >> _______________________________________________ >> python-control-discuss mailing list >> pyt...@li... >> https://lists.sourceforge.net/lists/listinfo/python-control-discuss >> > > ------------------------------------------------------------------------------ > The Next 800 Companies to Lead America's Growth: New Video Whitepaper > David G. Thomson, author of the best-selling book "Blueprint to a > Billion" shares his insights and actions to help propel your > business during the next growth cycle. Listen Now! > http://p.sf.net/sfu/SAP-dev2dev > _______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss |
|
From: Ryan K. <rk...@si...> - 2010-11-06 15:22:17
|
FYI, I saw some mention of lsim2 in a recent email. My code includes a step response function that first tries lsim then uses lsim2 if there is an exception. The exceptions usually are caused by pure integrators. I think that defaulting to lsim2 will make the simulation slower. I have fixed a few bugs in my code. There may even be a bug fix in the digital conversion code that hasn't made it into my main branch. The most recent versions of my code are on github and on my course website for this semester: https://github.com/ryanGT/controls http://www.cs.siue.edu/~rkrauss/mechatronics/2010/python_controls/python.html -- Ryan Krauss, Ph.D. Assistant Professor Mechanical Engineering Southern Illinois University Edwardsville On Sat, Nov 6, 2010 at 8:45 AM, Richard Murray <mu...@cd...> wrote: > A couple of updates on status of the the python-control package, related to the commit earlier today: > > * Ryan Krauss has sent me his control.py package that he uses for his class at Southern Illinois University Edwardsville. He has implemented functions for continuous to discrete time conversion and root locus plus. Ryan's code is now in trunk/external/control.py and Ryan has graciously given us permission to incorporate functions into the python-control package. > > * Roberto Bucher has sent me his yottalab.py package that has a number of functions for state space computations in continuous and discrete time. This code is in trunk/external/yottalab.py. Roberto has also give us permission to incoporate functions into the python-control library. > > * I've added a pade() command developed by Sawyer Fuller at Caltech, using the description in Golub and Van Loan. This code is in trunk/src/delay.py. > > * A group of students at Princeton (Lauren Padilla, Kevin Chen and Steve Brunton) are working on a class project this term and will be helping to implement some additional functions in the library. I've encouraged them to use this list to send out thoughts and ideas about their plans. > > * In integrating some of the functions from Ryan and Roberto, I've made various changes to things and fixed a couple of bugs. See trunk/ChangeLog for details. > > That's it for now. Comments and feedback is always welcome. > > -richard > > > > ------------------------------------------------------------------------------ > The Next 800 Companies to Lead America's Growth: New Video Whitepaper > David G. Thomson, author of the best-selling book "Blueprint to a > Billion" shares his insights and actions to help propel your > business during the next growth cycle. Listen Now! > http://p.sf.net/sfu/SAP-dev2dev > _______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss > |
|
From: Richard M. <mu...@cd...> - 2010-11-06 14:18:01
|
A couple of updates on status of the the python-control package, related to the commit earlier today: * Ryan Krauss has sent me his control.py package that he uses for his class at Southern Illinois University Edwardsville. He has implemented functions for continuous to discrete time conversion and root locus plus. Ryan's code is now in trunk/external/control.py and Ryan has graciously given us permission to incorporate functions into the python-control package. * Roberto Bucher has sent me his yottalab.py package that has a number of functions for state space computations in continuous and discrete time. This code is in trunk/external/yottalab.py. Roberto has also give us permission to incoporate functions into the python-control library. * I've added a pade() command developed by Sawyer Fuller at Caltech, using the description in Golub and Van Loan. This code is in trunk/src/delay.py. * A group of students at Princeton (Lauren Padilla, Kevin Chen and Steve Brunton) are working on a class project this term and will be helping to implement some additional functions in the library. I've encouraged them to use this list to send out thoughts and ideas about their plans. * In integrating some of the functions from Ryan and Roberto, I've made various changes to things and fixed a couple of bugs. See trunk/ChangeLog for details. That's it for now. Comments and feedback is always welcome. -richard |
|
From: <mur...@us...> - 2010-11-06 13:04:03
|
Revision: 29
http://python-control.svn.sourceforge.net/python-control/?rev=29&view=rev
Author: murrayrm
Date: 2010-11-06 13:03:55 +0000 (Sat, 06 Nov 2010)
Log Message:
-----------
2010-11-05 Richard Murray <murray@sumatra.local>
* external/yottalab.py: New file containing Roberto Bucher's control
library functions. OK to start pulling these into the main library,
with attribution, but note that they use modifications of the
default library => some rewrites will be needed.
2010-09-11 Richard Murray <murray@sumatra.local>
* src/matlab.py (step): Added local step response function that uses
lsim2() instead of signal.step (which can't handle integrators).
This function may not be needed when new scipy step2() function is
available.
(impulse): Added local impulse response function that sets the
initial condition based on the input matrix and then uses the
lsim2() function to compute the response.
* examples/test-response.py: Added test script for making sure that
time repsonse functions are working as desired
* src/matlab.py (lsim): Added local version of lsim that calls
signal.lsim2 (actual ODE integrator)
Modified Paths:
--------------
trunk/ChangeLog
trunk/Pending
trunk/src/__init__.py
trunk/src/matlab.py
trunk/src/rlocus.py
trunk/src/statefbk.py
Added Paths:
-----------
trunk/external/yottalab.py
trunk/src/delay.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2010-06-18 05:10:26 UTC (rev 28)
+++ trunk/ChangeLog 2010-11-06 13:03:55 UTC (rev 29)
@@ -1,3 +1,45 @@
+2010-11-05 Richard Murray <murray@sumatra.local>
+
+ * external/yottalab.py: New file containing Roberto Bucher's control
+ library functions. OK to start pulling these into the main library,
+ with attribution, but note that they use modifications of the
+ default library => some rewrites will be needed.
+
+2010-09-11 Richard Murray <murray@sumatra.local>
+
+ * src/matlab.py (step): Added local step response function that uses
+ lsim2() instead of signal.step (which can't handle integrators).
+ This function may not be needed when new scipy step2() function is
+ available.
+ (impulse): Added local impulse response function that sets the
+ initial condition based on the input matrix and then uses the
+ lsim2() function to compute the response.
+
+ * examples/test-response.py: Added test script for making sure that
+ time repsonse functions are working as desired
+
+ * src/matlab.py (lsim): Added local version of lsim that calls
+ signal.lsim2 (actual ODE integrator)
+
+2010-09-06 Richard Murray <murray@sumatra.local>
+
+ * src/statefbk.py (ctrb): new function for testing controllability
+ * src/statefbk.py (obsv): new function for testing observabiilty
+
+2010-09-02 Richard Murray <murray@sumatra.local>
+
+ * src/statefbk.py (place): Use np.size() instead of len() for
+ finding length of placed_eigs for better compatability with
+ different python versions [courtesy of Roberto Bucher]
+
+ * src/delay.py (pade): New file for delay-based computations +
+ initial implementation of pade() [courtesy Sawyer Fuller]
+
+2010-06-17 Richard Murray <murray@sumatra.local>
+
+ * src/rlocus.py: changed num, den to nump, denp for clarity
+ * src/rlocus.py: new file with Ryan Krauss's root locus code
+
2010-06-06 Richard Murray <murray@sumatra.local>
* examples/pvtol-lqr.py: Added example to test out LQR routines
Modified: trunk/Pending
===================================================================
--- trunk/Pending 2010-06-18 05:10:26 UTC (rev 28)
+++ trunk/Pending 2010-11-06 13:03:55 UTC (rev 29)
@@ -1,8 +1,42 @@
List of Pending changes for control-python
RMM, 5 Sep 09
+This file contains brief notes on features that need to be added to
+the python control library. Mainly intended to keep track of "bigger
+picture" things that need to be done.
+
+--> See src/matlab.py for a list of MATLAB functions that eventually need
+ to be implemented.
+
+OPEN BUGS
+ * step() doesn't handle systems with a pole at the origin (use lsim2)
+
+Transfer code from Roberto Bucher's yottalab to python-control
+ acker - pole placement using Ackermann method
+ c2d - contimous to discrete time conversion
+ full_obs - full order observer
+ red_obs - reduced order observer
+ comp_form - state feedback controller+observer in compact form
+ comp_form_i - state feedback controller+observer+integ in compact form
+ dsimul - simulate discrete time systems
+ dstep - step response (plot) of discrete time systems
+ dimpulse - imoulse response (plot) of discrete time systems
+ bb_step - step response (plot) of continous time systems
+ sysctr - system+controller+observer+feedback
+ care - Solve Riccati equation for contimous time systems
+ dare - Solve Riccati equation for discrete time systems
+ dlqr - discrete linear quadratic regulator
+ minreal - minimal state space representation
+
+Transfer code from Ryan Krauss's control.py to python-control
+ * phase margin computations (as part of margin command)
+ * step reponse
+ * c2d, c2d_tustin (compare to Bucher version first)
+
Examples and test cases
+ * Put together unit tests for all functions (after deciding on framework)
* Figure out how to import 'figure' command properly (version issue?)
+ * Figure out source of BadCoefficients warning messages (pvtol-lqr and others)
TransferFunction class fixes
* evalfr is not working (num, den stored as ndarrays, not poly1ds)
@@ -16,12 +50,11 @@
* Implement pzmap for state space systems
LTI updates
- * Try to get control.matlab.step semantics to be closer to MATLAB
+ * Implement control.matlab.step (with semantics similar to MATLAB)
Basic functions to be added
* margin - compute gain and phase margin (no plot)
* lqr - compute optimal feedback gains (use SLICOT SB02ND.f)
- * lqe - compute optimal feedback gains (use SLICOT SB02ND.f)
* lyap - solve Lyapunov equation (use SLICOT SB03MD.f)
* See http://www.slicot.org/shared/libindex.html for list of functions
Added: trunk/external/yottalab.py
===================================================================
--- trunk/external/yottalab.py (rev 0)
+++ trunk/external/yottalab.py 2010-11-06 13:03:55 UTC (rev 29)
@@ -0,0 +1,652 @@
+"""
+This is a procedural interface to the yttalab library
+
+ro...@su...
+
+The following commands are provided:
+
+Design and plot commands
+ acker - pole placement using Ackermann method
+ c2d - contimous to discrete time conversion
+ full_obs - full order observer
+ red_obs - reduced order observer
+ comp_form - state feedback controller+observer in compact form
+ comp_form_i - state feedback controller+observer+integ in compact form
+ dsimul - simulate discrete time systems
+ dstep - step response (plot) of discrete time systems
+ dimpulse - imoulse response (plot) of discrete time systems
+ bb_step - step response (plot) of continous time systems
+ sysctr - system+controller+observer+feedback
+ care - Solve Riccati equation for contimous time systems
+ dare - Solve Riccati equation for discrete time systems
+ dlqr - discrete linear quadratic regulator
+ minreal - minimal state space representation
+
+"""
+from matplotlib.pylab import *
+from control.matlab import *
+from numpy import hstack,vstack,pi
+from scipy import zeros,ones,eye,mat,shape,size,size, \
+ arange,real,poly,array,diag
+from scipy.linalg import det,inv,expm,eig,eigvals
+import numpy as np
+import scipy as sp
+from slycot import sb02od, tb03ad
+
+def acker(A,B,poles):
+ """Pole placemenmt using Ackermann method
+
+ Call:
+ k=acker(A,B,poles)
+
+ Parameters
+ ----------
+ A, B : State and input matrix of the system
+ poles: desired poles
+
+ Returns
+ -------
+ k: matrix
+ State feedback gains
+
+ """
+ a=mat(A)
+ b=mat(B)
+ p=real(poly(poles))
+ ct=ctrb(A,B)
+ if det(ct)==0:
+ k=0
+ print "Pole placement invalid"
+ else:
+ n=size(p)
+ pmat=p[n-1]*a**0
+ for i in arange(1,n):
+ pmat=pmat+p[n-i-1]*a**i
+ k=inv(ct)*pmat
+ k=k[-1][:]
+ return k
+
+def c2d(sys,Ts):
+ """Continous to discrete conversion with ZOH method
+
+ Call:
+ sysd=c2d(sys,Ts)
+
+ Parameters
+ ----------
+ sys : System in statespace or Tf form
+ Ts: Sampling Time
+
+ Returns
+ -------
+ sysd: ss or Tf system
+ Discrete system
+
+ """
+ flag = 0
+ if type(sys).__name__=='TransferFunction':
+ sys=tf2ss(sys)
+ flag=1
+ a=sys.A
+ b=sys.B
+ c=sys.C
+ d=sys.D
+ n=shape(a)[0]
+ nb=shape(b)[1]
+ ztmp=zeros((nb,n+nb))
+ tmp=hstack((a,b))
+ tmp=vstack((tmp,ztmp))
+ tmp=expm(tmp*Ts)
+ a=tmp[0:n,0:n]
+ b=tmp[0:n,n:n+nb]
+ sysd=ss(a,b,c,d,Ts)
+ if flag==1:
+ sysd=ss2tf(sysd)
+ return sysd
+
+def full_obs(sys,poles):
+ """Full order observer of the system sys
+
+ Call:
+ obs=full_obs(sys,poles)
+
+ Parameters
+ ----------
+ sys : System in State Space form
+ poles: desired observer poles
+
+ Returns
+ -------
+ obs: ss
+ Observer
+
+ """
+ if type(sys).__name__=='TransferFunction':
+ "System must be in state space form"
+ return
+ a=mat(sys.A)
+ b=mat(sys.B)
+ c=mat(sys.C)
+ d=mat(sys.D)
+ poles=mat(poles)
+ L=place(a.T,c.T,poles)
+ L=mat(L).T
+ Ao=a-L*c
+ Bo=hstack((b-L*d,L))
+ n=shape(Ao)
+ m=shape(Bo)
+ Co=eye(n[0],n[1])
+ Do=zeros((n[0],m[1]))
+ obs=ss(Ao,Bo,Co,Do,sys.Tsamp)
+ return obs
+
+def red_obs(sys,T,poles):
+ """Reduced order observer of the system sys
+
+ Call:
+ obs=red_obs(sys,T,poles)
+
+ Parameters
+ ----------
+ sys : System in State Space form
+ T: Complement matrix
+ poles: desired observer poles
+
+ Returns
+ -------
+ obs: ss
+ Reduced order Observer
+
+ """
+ if type(sys).__name__=='TransferFunction':
+ "System must be in state space form"
+ return
+ a=mat(sys.A)
+ b=mat(sys.B)
+ c=mat(sys.C)
+ d=mat(sys.D)
+ T=mat(T)
+ P=mat(vstack((c,T)))
+ poles=mat(poles)
+ invP=inv(P)
+ AA=P*a*invP
+ ny=shape(c)[0]
+ nx=shape(a)[0]
+ nu=shape(b)[1]
+
+ A11=AA[0:ny,0:ny]
+ A12=AA[0:ny,ny:nx]
+ A21=AA[ny:nx,0:ny]
+ A22=AA[ny:nx,ny:nx]
+
+ L1=place(A22.T,A12.T,poles)
+ L1=mat(L1).T
+
+ nn=nx-ny
+
+ tmp1=mat(hstack((-L1,eye(nn,nn))))
+ tmp2=mat(vstack((zeros((ny,nn)),eye(nn,nn))))
+ Ar=tmp1*P*a*invP*tmp2
+
+ tmp3=vstack((eye(ny,ny),L1))
+ tmp3=mat(hstack((P*b,P*a*invP*tmp3)))
+ tmp4=hstack((eye(nu,nu),zeros((nu,ny))))
+ tmp5=hstack((-d,eye(ny,ny)))
+ tmp4=mat(vstack((tmp4,tmp5)))
+
+ Br=tmp1*tmp3*tmp4
+
+ Cr=invP*tmp2
+
+ tmp5=hstack((zeros((ny,nu)),eye(ny,ny)))
+ tmp6=hstack((zeros((nn,nu)),L1))
+ tmp5=mat(vstack((tmp5,tmp6)))
+ Dr=invP*tmp5*tmp4
+
+ obs=ss(Ar,Br,Cr,Dr,sys.Tsamp)
+ return obs
+
+def comp_form(sys,obs,K):
+ """Compact form Conroller+Observer
+
+ Call:
+ contr=comp_form(sys,obs,K)
+
+ Parameters
+ ----------
+ sys : System in State Space form
+ obs : Observer in State Space form
+ K: State feedback gains
+
+ Returns
+ -------
+ contr: ss
+ Controller
+
+ """
+ nx=shape(sys.A)[0]
+ ny=shape(sys.C)[0]
+ nu=shape(sys.B)[1]
+ no=shape(obs.A)[0]
+
+ Bu=mat(obs.B[:,0:nu])
+ By=mat(obs.B[:,nu:])
+ Du=mat(obs.D[:,0:nu])
+ Dy=mat(obs.D[:,nu:])
+
+ X=inv(eye(nu,nu)+K*Du)
+
+ Ac = mat(obs.A)-Bu*X*K*mat(obs.C);
+ Bc = hstack((Bu*X,By-Bu*X*K*Dy))
+ Cc = -X*K*mat(obs.C);
+ Dc = hstack((X,-X*K*Dy))
+ contr = ss(Ac,Bc,Cc,Dc,sys.Tsamp)
+ return contr
+
+def comp_form_i(sys,obs,K,Ts,Cy=[[1]]):
+ """Compact form Conroller+Observer+Integral part
+ Only for discrete systems!!!
+
+ Call:
+ contr=comp_form_i(sys,obs,K,Ts[,Cy])
+
+ Parameters
+ ----------
+ sys : System in State Space form
+ obs : Observer in State Space form
+ K: State feedback gains
+ Ts: Sampling time
+ Cy: feedback matric to choose the output for integral part
+
+ Returns
+ -------
+ contr: ss
+ Controller
+
+ """
+ if sys.Tsamp==0.0:
+ print "contr_form_i works only with discrete systems!"
+ return
+
+ ny=shape(sys.C)[0]
+ nu=shape(sys.B)[1]
+ nx=shape(sys.A)[0]
+ no=shape(obs.A)[0]
+ ni=shape(Cy)[0]
+
+ B_obsu = mat(obs.B[:,0:nu])
+ B_obsy = mat(obs.B[:,nu:nu+ny])
+ D_obsu = mat(obs.D[:,0:nu])
+ D_obsy = mat(obs.D[:,nu:nu+ny])
+
+ k=mat(K)
+ nk=shape(k)[1]
+ Ke=k[:,nk-ni:]
+ K=k[:,0:nk-ni]
+ X = inv(eye(nu,nu)+K*D_obsu);
+
+ a=mat(obs.A)
+ c=mat(obs.C)
+ Cy=mat(Cy)
+
+ tmp1=hstack((a-B_obsu*X*K*c,-B_obsu*X*Ke))
+
+ tmp2=hstack((zeros((ni,no)),eye(ni,ni)))
+ A_ctr=vstack((tmp1,tmp2))
+
+ tmp1=hstack((zeros((no,ni)),-B_obsu*X*K*D_obsy+B_obsy))
+ tmp2=hstack((eye(ni,ni)*Ts,-Cy*Ts))
+ B_ctr=vstack((tmp1,tmp2))
+
+ C_ctr=hstack((-X*K*c,-X*Ke))
+ D_ctr=hstack((zeros((nu,ni)),-X*K*D_obsy))
+
+ contr=ss(A_ctr,B_ctr,C_ctr,D_ctr,sys.Tsamp)
+ return contr
+
+def dsimul(sys,u):
+ """Simulate the discrete system sys
+ Only for discrete systems!!!
+
+ Call:
+ y=dsimul(sys,u)
+
+ Parameters
+ ----------
+ sys : Discrete System in State Space form
+ u : input vector
+ Returns
+ -------
+ y: ndarray
+ Simulation results
+
+ """
+ a=mat(sys.A)
+ b=mat(sys.B)
+ c=mat(sys.C)
+ d=mat(sys.D)
+ nx=shape(a)[0]
+ ns=shape(u)[1]
+ xk=zeros((nx,1))
+ for i in arange(0,ns):
+ uk=u[:,i]
+ xk_1=a*xk+b*uk
+ yk=c*xk+d*uk
+ xk=xk_1
+ if i==0:
+ y=yk
+ else:
+ y=hstack((y,yk))
+ y=array(y).T
+ return y
+
+def dstep(sys,Tf=10.0):
+ """Plot the step response of the discrete system sys
+ Only for discrete systems!!!
+
+ Call:
+ y=dstep(sys, [,Tf=final time]))
+
+ Parameters
+ ----------
+ sys : Discrete System in State Space form
+ Tf : Final simulation time
+
+ Returns
+ -------
+ Nothing
+
+ """
+ Ts=sys.Tsamp
+ if Ts==0.0:
+ "Only discrete systems allowed!"
+ return
+
+ ns=int(Tf/Ts+1)
+ u=ones((1,ns))
+ y=dsimul(sys,u)
+ T=arange(0,Tf+Ts/2,Ts)
+ plot(T,y)
+ grid()
+ show()
+
+def dimpulse(sys,Tf=10.0):
+ """Plot the impulse response of the discrete system sys
+ Only for discrete systems!!!
+
+ Call:
+ y=dimpulse(sys,[,Tf=final time]))
+
+ Parameters
+ ----------
+ sys : Discrete System in State Space form
+ Tf : Final simulation time
+
+ Returns
+ -------
+ Nothing
+
+ """
+ Ts=sys.Tsamp
+ if Ts==0.0:
+ "Only discrete systems allowed!"
+ return
+
+ ns=int(Tf/Ts+1)
+ u=zeros((1,ns))
+ u[0,0]=1/Ts
+ y=dsimul(sys,u)
+ T=arange(0,Tf+Ts/2,Ts)
+ plot(T,y)
+ grid()
+ show()
+
+# Step response (plot)
+def bb_step(sys,X0=None,Tf=None,Ts=0.001):
+ """Plot the step response of the continous system sys
+
+ Call:
+ y=bb_step(sys [,Tf=final time] [,Ts=time step])
+
+ Parameters
+ ----------
+ sys : Continous System in State Space form
+ X0: Initial state vector (not used yet)
+ Ts : sympling time
+ Tf : Final simulation time
+
+ Returns
+ -------
+ Nothing
+
+ """
+ if Tf==None:
+ vals = eigvals(sys.A)
+ r = min(abs(real(vals)))
+ if r == 0.0:
+ r = 1.0
+ Tf = 7.0 / r
+ sysd=c2d(sys,Ts)
+ dstep(sysd,Tf=Tf)
+
+def sysctr(sys,contr):
+ """Build the discrete system controller+plant+output feedback
+
+ Call:
+ syscontr=sysctr(sys,contr)
+
+ Parameters
+ ----------
+ sys : Continous System in State Space form
+ contr: Controller (with observer if required)
+
+ Returns
+ -------
+ sysc: ss system
+ The system with reference as input and outputs of plants
+ as output
+
+ """
+ if contr.Tsamp!=sys.Tsamp:
+ print "Systems with different sampling time!!!"
+ return
+ sysf=contr*sys
+
+ nu=shape(sysf.B)[1]
+ b1=mat(sysf.B[:,0])
+ b2=mat(sysf.B[:,1:nu])
+ d1=mat(sysf.D[:,0])
+ d2=mat(sysf.D[:,1:nu])
+
+ n2=shape(d2)[0]
+
+ Id=mat(eye(n2,n2))
+ X=inv(Id-d2)
+
+ Af=mat(sysf.A)+b2*X*mat(sysf.C)
+ Bf=b1+b2*X*d1
+ Cf=X*mat(sysf.C)
+ Df=X*d1
+
+ sysc=ss(Af,Bf,Cf,Df,sys.Tsamp)
+ return sysc
+
+def care(A,B,Q,R):
+ """Solve Riccati equation for discrete time systems
+
+ Usage
+ =====
+ [K, S, E] = care(A, B, Q, R)
+
+ Inputs
+ ------
+ A, B: 2-d arrays with dynamics and input matrices
+ sys: linear I/O system
+ Q, R: 2-d array with state and input weight matrices
+
+ Outputs
+ -------
+ X: solution of the Riccati eq.
+ """
+
+ # Check dimensions for consistency
+ nstates = B.shape[0];
+ ninputs = B.shape[1];
+ if (A.shape[0] != nstates or A.shape[1] != nstates):
+ raise ControlDimension("inconsistent system dimensions")
+
+ elif (Q.shape[0] != nstates or Q.shape[1] != nstates or
+ R.shape[0] != ninputs or R.shape[1] != ninputs) :
+ raise ControlDimension("incorrect weighting matrix dimensions")
+
+ rcond,X,w,S,T = \
+ sb02od(nstates, ninputs, A, B, Q, R, 'C');
+
+ return X
+
+
+def dare(A,B,Q,R):
+ """Solve Riccati equation for discrete time systems
+
+ Usage
+ =====
+ [K, S, E] = care(A, B, Q, R)
+
+ Inputs
+ ------
+ A, B: 2-d arrays with dynamics and input matrices
+ sys: linear I/O system
+ Q, R: 2-d array with state and input weight matrices
+
+ Outputs
+ -------
+ X: solution of the Riccati eq.
+ """
+
+ # Check dimensions for consistency
+ nstates = B.shape[0];
+ ninputs = B.shape[1];
+ if (A.shape[0] != nstates or A.shape[1] != nstates):
+ raise ControlDimension("inconsistent system dimensions")
+
+ elif (Q.shape[0] != nstates or Q.shape[1] != nstates or
+ R.shape[0] != ninputs or R.shape[1] != ninputs) :
+ raise ControlDimension("incorrect weighting matrix dimensions")
+
+ rcond,X,w,S,T = \
+ sb02od(nstates, ninputs, A, B, Q, R, 'D');
+
+ return X
+
+
+def dlqr(*args, **keywords):
+ """Linear quadratic regulator design for discrete systems
+
+ Usage
+ =====
+ [K, S, E] = dlqr(A, B, Q, R, [N])
+ [K, S, E] = dlqr(sys, Q, R, [N])
+
+ The dlqr() function computes the optimal state feedback controller
+ that minimizes the quadratic cost
+
+ J = \sum_0^\infty x' Q x + u' R u + 2 x' N u
+
+ Inputs
+ ------
+ A, B: 2-d arrays with dynamics and input matrices
+ sys: linear I/O system
+ Q, R: 2-d array with state and input weight matrices
+ N: optional 2-d array with cross weight matrix
+
+ Outputs
+ -------
+ K: 2-d array with state feedback gains
+ S: 2-d array with solution to Riccati equation
+ E: 1-d array with eigenvalues of the closed loop system
+ """
+
+ #
+ # Process the arguments and figure out what inputs we received
+ #
+
+ # Get the system description
+ if (len(args) < 3):
+ raise ControlArgument("not enough input arguments")
+
+ elif (ctrlutil.issys(args[0])):
+ # We were passed a system as the first argument; extract A and B
+ A = array(args[0].A, ndmin=2, dtype=float);
+ B = array(args[0].B, ndmin=2, dtype=float);
+ index = 1;
+ if args[0].Tsamp==0.0:
+ print "dlqr works only for discrete systems!"
+ return
+ else:
+ # Arguments should be A and B matrices
+ A = array(args[0], ndmin=2, dtype=float);
+ B = array(args[1], ndmin=2, dtype=float);
+ index = 2;
+
+ # Get the weighting matrices (converting to matrices, if needed)
+ Q = array(args[index], ndmin=2, dtype=float);
+ R = array(args[index+1], ndmin=2, dtype=float);
+ if (len(args) > index + 2):
+ N = array(args[index+2], ndmin=2, dtype=float);
+ else:
+ N = zeros((Q.shape[0], R.shape[1]));
+
+ # Check dimensions for consistency
+ nstates = B.shape[0];
+ ninputs = B.shape[1];
+ if (A.shape[0] != nstates or A.shape[1] != nstates):
+ raise ControlDimension("inconsistent system dimensions")
+
+ elif (Q.shape[0] != nstates or Q.shape[1] != nstates or
+ R.shape[0] != ninputs or R.shape[1] != ninputs or
+ N.shape[0] != nstates or N.shape[1] != ninputs):
+ raise ControlDimension("incorrect weighting matrix dimensions")
+
+ #Solve the riccati equation
+ X = dare(A,B,Q,R)
+
+ # Now compute the return value
+ Phi=mat(A)
+ H=mat(B)
+ K=inv(H.T*X*H+R)*(H.T*X*Phi)
+ L=eig(Phi-H*K)
+ return K,X,L
+
+def minreal(sys):
+ """Minimal representation for state space systems
+
+ Usage
+ =====
+ [sysmin]=minreal[sys]
+
+ Inputs
+ ------
+
+ sys: system in ss or tf form
+
+ Outputs
+ -------
+ sysfin: system in state space form
+ """
+ a=mat(sys.A)
+ b=mat(sys.B)
+ c=mat(sys.C)
+ d=mat(sys.D)
+ nx=shape(a)[0]
+ ni=shape(b)[1]
+ no=shape(c)[0]
+
+ out=tb03ad(nx,no,ni,a,b,c,d,'R')
+
+ nr=out[3]
+ A=out[0][:nr,:nr]
+ B=out[1][:nr,:ni]
+ C=out[2][:no,:nr]
+ sysf=ss(A,B,C,sys.D,sys.Tsamp)
+ return sysf
+
Modified: trunk/src/__init__.py
===================================================================
--- trunk/src/__init__.py 2010-06-18 05:10:26 UTC (rev 28)
+++ trunk/src/__init__.py 2010-11-06 13:03:55 UTC (rev 29)
@@ -61,3 +61,4 @@
from freqplot import *
from bdalg import *
from statefbk import *
+from delay import *
Added: trunk/src/delay.py
===================================================================
--- trunk/src/delay.py (rev 0)
+++ trunk/src/delay.py 2010-11-06 13:03:55 UTC (rev 29)
@@ -0,0 +1,68 @@
+# delay.py - functions involving time delays
+#
+# Author: Sawyer Fuller
+# Date: 26 Aug 2010
+#
+# This file contains functions for implementing time delays (currently
+# only the pade() function).
+#
+# Copyright (c) 2010 by California Institute of Technology
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions
+# are met:
+#
+# 1. Redistributions of source code must retain the above copyright
+# notice, this list of conditions and the following disclaimer.
+#
+# 2. Redistributions in binary form must reproduce the above copyright
+# notice, this list of conditions and the following disclaimer in the
+# documentation and/or other materials provided with the distribution.
+#
+# 3. Neither the name of the California Institute of Technology nor
+# the names of its contributors may be used to endorse or promote
+# products derived from this software without specific prior
+# written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CALTECH
+# OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
+# USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+# OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+# SUCH DAMAGE.
+#
+# $Id: pade.py 17 2010-05-29 23:50:52Z murrayrm $
+
+def pade(T, n=1):
+ """
+ Return the numerator and denominator coefficients of the Pade approximation.
+
+ Inputs:
+ T --> time delay
+ n --> order of approximation
+
+ Outputs:
+ num, den --> arrays in descending powers of s.
+
+ Based on an algorithm in Golub and van Loan, "Matrix Computation" 3rd.
+ Ed. pp. 572-574.
+ """
+ num = [0. for i in range(n+1)]
+ num[-1] = 1.
+ den = [0. for i in range(n+1)]
+ den[-1] = 1.
+ c = 1.
+ for k in range(1, n+1):
+ c = T * c * (n - k + 1)/(2 * n - k + 1)/k
+ num[n - k] = c * (-1)**k
+ den[n - k] = c
+ num = [coeff/den[0] for coeff in num]
+ den = [coeff/den[0] for coeff in den]
+ return num, den
Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py 2010-06-18 05:10:26 UTC (rev 28)
+++ trunk/src/matlab.py 2010-11-06 13:03:55 UTC (rev 29)
@@ -55,7 +55,6 @@
# Import MATLAB-like functions that are defined in other packages
from scipy.signal import zpk2ss, ss2zpk, tf2zpk, zpk2tf
-from scipy.signal import lsim, impulse, step
from scipy import linspace, logspace
# Control system library
@@ -70,7 +69,8 @@
from freqplot import nyquist, gangof4
from bdalg import series, parallel, negate, feedback
from pzmap import pzmap
-from statefbk import place, lqr
+from statefbk import ctrb, obsv, place, lqr
+from delay import pade
__doc__ = """
The control.matlab module defines functions that are roughly the
@@ -191,8 +191,8 @@
drss - random stable discrete-time state-space models
ss2ss - state coordinate transformation
canon - canonical forms of state-space models
- ctrb - controllability matrix
- obsv - observability matrix
+* ctrb - controllability matrix
+* obsv - observability matrix
gram - controllability and observability gramians
ss/prescale - optimal scaling of state-space models.
balreal - gramian-based input/output balancing
@@ -214,7 +214,7 @@
lti/hasdelay - true for models with time delays
lti/totaldelay - total delay between each input/output pair
lti/delay2z - replace delays by poles at z=0 or FRD phase shift
- pade - pade approximation of time delays
+* pade - pade approximation of time delays
Model dimensions and characteristics
class - model type ('tf', 'zpk', 'ss', or 'frd')
@@ -370,3 +370,89 @@
# Call the bode command
return freqplot.bode(syslist, omega, **keywords)
+
+#
+# Modifications to scipy.signal functions
+#
+
+# Redefine lsim to use lsim2
+def lsim(*args, **keywords):
+ """Simulate the output of a linear system
+
+ Usage
+ =====
+ (T, yout, xout) = lsim(sys, u, T, X0)
+
+ Inputs:
+ sys LTI system
+ u input array giving input at each time T
+ T time steps at which the input is defined
+ X0 initial condition (optional, default = 0)
+
+ Outputs:
+ T time values of the output
+ yout response of the system
+ xout time evolution of the state vector
+ """
+ return sp.signal.lsim2(*args, **keywords)
+
+#! Redefine step to use lsim2
+#! Not yet implemented
+def step(*args, **keywords):
+ """Step response of a linear system
+
+ Usage
+ =====
+ (T, yout) = step(sys, T, X0)
+
+ Inputs:
+ sys LTI system
+ T time steps (optional; autocomputed if not gien)
+ X0 initial condition (optional, default = 0)
+
+ Outputs:
+ T time values of the output
+ yout response of the system
+ """
+ return sp.signal.step(*args, **keywords)
+
+# Redefine initial to use lsim2
+#! Not yet implemented (uses step for now)
+def initial(*args, **keywords):
+ """Initial condition response of a linear system
+
+ Usage
+ =====
+ (T, yout) = initial(sys, T, X0)
+
+ Inputs:
+ sys LTI system
+ T time steps (optional; autocomputed if not gien)
+ X0 initial condition (optional, default = 0)
+
+ Outputs:
+ T time values of the output
+ yout response of the system
+ """
+ return sp.signal.initial(*args, **keywords)
+
+# Redefine impulse to use initial()
+#! Not yet implemented (uses impulse for now)
+def impulse(*args, **keywords):
+ """Step response of a linear system
+
+ Usage
+ =====
+ (T, yout) = impulse(sys, T, X0)
+
+ Inputs:
+ sys LTI system
+ T time steps (optional; autocomputed if not gien)
+ X0 initial condition (optional, default = 0)
+
+ Outputs:
+ T time values of the output
+ yout response of the system
+ """
+ return sp.signal.impulse(*args, **keywords)
+
Modified: trunk/src/rlocus.py
===================================================================
--- trunk/src/rlocus.py 2010-06-18 05:10:26 UTC (rev 28)
+++ trunk/src/rlocus.py 2010-11-06 13:03:55 UTC (rev 29)
@@ -53,7 +53,7 @@
of kvect."""
# Convert numerator and denominator to polynomials if they aren't
- (num, den) = _systopoly1d(sys);
+ (nump, denp) = _systopoly1d(sys);
# Set up the figure
if fig is None:
@@ -68,11 +68,11 @@
mymat = _RLSortRoots(sys, mymat)
# plot open loop poles
- poles = array(den.r)
+ poles = array(denp.r)
ax.plot(real(poles), imag(poles), 'x')
# plot open loop zeros
- zeros = array(num.r)
+ zeros = array(nump.r)
if zeros.any():
ax.plot(real(zeros), imag(zeros), 'o')
@@ -94,23 +94,23 @@
"""Extract numerator and denominator polynomails for a system"""
# Start by extracting the numerator and denominator from system object
- num = sys.num; den = sys.den;
+ nump = sys.num; denp = sys.den;
# Check to see if num, den are already polynomials; otherwise convert
- if (not isinstance(num, poly1d)): num = poly1d(num)
- if (not isinstance(den, poly1d)): den = poly1d(den)
+ if (not isinstance(nump, poly1d)): nump = poly1d(nump)
+ if (not isinstance(denp, poly1d)): denp = poly1d(denp)
- return (num, den)
+ return (nump, denp)
def _RLFindRoots(sys, kvect):
"""Find the roots for the root locus."""
# Convert numerator and denominator to polynomials if they aren't
- (num, den) = _systopoly1d(sys);
+ (nump, denp) = _systopoly1d(sys);
roots = []
for k in kvect:
- curpoly = den + k * num
+ curpoly = denp + k * nump
curroots = curpoly.r
curroots.sort()
roots.append(curroots)
Modified: trunk/src/statefbk.py
===================================================================
--- trunk/src/statefbk.py 2010-06-18 05:10:26 UTC (rev 28)
+++ trunk/src/statefbk.py 2010-11-06 13:03:55 UTC (rev 29)
@@ -1,6 +1,6 @@
# statefbk.py - tools for state feedback control
#
-# Author: Richard M. Murray
+# Author: Richard M. Murray, Roberto Bucher
# Date: 31 May 2010
#
# This file contains routines for designing state space controllers
@@ -92,7 +92,7 @@
# Call SLICOT routine to place the eigenvalues
A_z,w,nfp,nap,nup,F,Z = \
- sb01bd(B_mat.shape[0], B_mat.shape[1], len(placed_eigs), alpha,
+ sb01bd(B_mat.shape[0], B_mat.shape[1], np.size(placed_eigs), alpha,
A_mat, B_mat, placed_eigs, 'C');
# Return the gain matrix, with MATLAB gain convention
@@ -184,3 +184,57 @@
E = w[0:nstates];
return K, S, E
+
+def ctrb(A,B):
+ """Controllabilty matrix
+
+ Usage
+ =====
+ Wc = ctrb(A, B)
+
+ Inputs
+ ------
+ A, B: Dynamics and input matrix of the system
+
+ Outputs
+ -------
+ Wc: Controllability matrix
+ """
+
+ # Convert input parameters to matrices (if they aren't already)
+ amat = np.mat(A)
+ bmat = np.mat(B)
+ n = np.shape(amat)[0]
+
+ # Construct the controllability matrix
+ ctrb = bmat
+ for i in range(1, n):
+ ctrb = np.vstack((ctrb, amat**i*bmat))
+ return ctrb
+
+def obsv(A, C):
+ """Observability matrix
+
+ Usage
+ =====
+ Wc = obsv(A, C)
+
+ Inputs
+ ------
+ A, C: Dynamics and output matrix of the system
+
+ Outputs
+ -------
+ Wc: Observability matrix
+ """
+
+ # Convert input parameters to matrices (if they aren't already)
+ amat = np.mat(A)
+ cmat = np.mat(C)
+ n = np.shape(amat)[0]
+
+ # Construct the controllability matrix
+ obsv = cmat
+ for i in range(1, n):
+ obsv = np.hstack((obsv, cmat*amat**i))
+ return obsv
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mur...@us...> - 2010-06-18 05:10:32
|
Revision: 28
http://python-control.svn.sourceforge.net/python-control/?rev=28&view=rev
Author: murrayrm
Date: 2010-06-18 05:10:26 +0000 (Fri, 18 Jun 2010)
Log Message:
-----------
set svn:keywords to Id to get tag working
Modified Paths:
--------------
trunk/external/controls.py
trunk/src/rlocus.py
Property Changed:
----------------
trunk/external/controls.py
trunk/src/rlocus.py
Modified: trunk/external/controls.py
===================================================================
--- trunk/external/controls.py 2010-06-18 05:06:29 UTC (rev 27)
+++ trunk/external/controls.py 2010-06-18 05:10:26 UTC (rev 28)
@@ -1,5 +1,5 @@
# controls.py - Ryan Krauss's control module
-# $Id: $
+# $Id$
"""This module is for analyzing linear, time-invariant dynamic systems
and feedback control systems using the Laplace transform. The heart
Property changes on: trunk/external/controls.py
___________________________________________________________________
Added: svn:keywords
+ Id
Modified: trunk/src/rlocus.py
===================================================================
--- trunk/src/rlocus.py 2010-06-18 05:06:29 UTC (rev 27)
+++ trunk/src/rlocus.py 2010-06-18 05:10:26 UTC (rev 28)
@@ -40,7 +40,7 @@
# or a controls.TransferFunction object.
# * Added some comments to make sure I understand the code
#
-# $Id: $
+# $Id$
# Packages used by this module
from scipy import *
Property changes on: trunk/src/rlocus.py
___________________________________________________________________
Added: svn:keywords
+ Id
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mur...@us...> - 2010-06-18 05:06:36
|
Revision: 27
http://python-control.svn.sourceforge.net/python-control/?rev=27&view=rev
Author: murrayrm
Date: 2010-06-18 05:06:29 +0000 (Fri, 18 Jun 2010)
Log Message:
-----------
Extracted Ryan's root locus code to see how things work. Main change
was converting num, den to poly1d's (assumed in Ryan's TF class).
Basic functionality seems to work (!).
Modified Paths:
--------------
trunk/external/controls.py
Added Paths:
-----------
trunk/src/rlocus.py
Modified: trunk/external/controls.py
===================================================================
--- trunk/external/controls.py 2010-06-17 21:29:23 UTC (rev 26)
+++ trunk/external/controls.py 2010-06-18 05:06:29 UTC (rev 27)
@@ -1,5 +1,5 @@
# controls.py - Ryan Krauss's control module
-# $Id$
+# $Id: $
"""This module is for analyzing linear, time-invariant dynamic systems
and feedback control systems using the Laplace transform. The heart
Added: trunk/src/rlocus.py
===================================================================
--- trunk/src/rlocus.py (rev 0)
+++ trunk/src/rlocus.py 2010-06-18 05:06:29 UTC (rev 27)
@@ -0,0 +1,140 @@
+# rlocus.py - code for computing a root locus plot
+# Code contributed by Ryan Krauss, 2010
+#
+# Copyright (c) 2010 by Ryan Krauss
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions
+# are met:
+#
+# 1. Redistributions of source code must retain the above copyright
+# notice, this list of conditions and the following disclaimer.
+#
+# 2. Redistributions in binary form must reproduce the above copyright
+# notice, this list of conditions and the following disclaimer in the
+# documentation and/or other materials provided with the distribution.
+#
+# 3. Neither the name of the California Institute of Technology nor
+# the names of its contributors may be used to endorse or promote
+# products derived from this software without specific prior
+# written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CALTECH
+# OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
+# USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+# OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+# SUCH DAMAGE.
+#
+# RMM, 17 June 2010: modified to be a standalone piece of code
+# * Added BSD copyright info to file (per Ryan)
+# * Added code to convert (num, den) to poly1d's if they aren't already.
+# This allows Ryan's code to run on a standard signal.ltisys object
+# or a controls.TransferFunction object.
+# * Added some comments to make sure I understand the code
+#
+# $Id: $
+
+# Packages used by this module
+from scipy import *
+
+# Main function: compute a root locus diagram
+def RootLocus(sys, kvect, fig=None, fignum=1, \
+ clear=True, xlim=None, ylim=None, plotstr='-'):
+ """Calculate the root locus by finding the roots of 1+k*TF(s)
+ where TF is self.num(s)/self.den(s) and each k is an element
+ of kvect."""
+
+ # Convert numerator and denominator to polynomials if they aren't
+ (num, den) = _systopoly1d(sys);
+
+ # Set up the figure
+ if fig is None:
+ import pylab
+ fig = pylab.figure(fignum)
+ if clear:
+ fig.clf()
+ ax = fig.add_subplot(111)
+
+ # Compute out the loci
+ mymat = _RLFindRoots(sys, kvect)
+ mymat = _RLSortRoots(sys, mymat)
+
+ # plot open loop poles
+ poles = array(den.r)
+ ax.plot(real(poles), imag(poles), 'x')
+
+ # plot open loop zeros
+ zeros = array(num.r)
+ if zeros.any():
+ ax.plot(real(zeros), imag(zeros), 'o')
+
+ # Now plot the loci
+ for col in mymat.T:
+ ax.plot(real(col), imag(col), plotstr)
+
+ # Set up plot axes and labels
+ if xlim:
+ ax.set_xlim(xlim)
+ if ylim:
+ ax.set_ylim(ylim)
+ ax.set_xlabel('Real')
+ ax.set_ylabel('Imaginary')
+ return mymat
+
+# Utility function to extract numerator and denominator polynomials
+def _systopoly1d(sys):
+ """Extract numerator and denominator polynomails for a system"""
+
+ # Start by extracting the numerator and denominator from system object
+ num = sys.num; den = sys.den;
+
+ # Check to see if num, den are already polynomials; otherwise convert
+ if (not isinstance(num, poly1d)): num = poly1d(num)
+ if (not isinstance(den, poly1d)): den = poly1d(den)
+
+ return (num, den)
+
+def _RLFindRoots(sys, kvect):
+ """Find the roots for the root locus."""
+
+ # Convert numerator and denominator to polynomials if they aren't
+ (num, den) = _systopoly1d(sys);
+
+ roots = []
+ for k in kvect:
+ curpoly = den + k * num
+ curroots = curpoly.r
+ curroots.sort()
+ roots.append(curroots)
+ mymat = row_stack(roots)
+ return mymat
+
+def _RLSortRoots(sys, mymat):
+ """Sort the roots from sys._RLFindRoots, so that the root
+ locus doesn't show weird pseudo-branches as roots jump from
+ one branch to another."""
+
+ sorted = zeros_like(mymat)
+ for n, row in enumerate(mymat):
+ if n==0:
+ sorted[n,:] = row
+ else:
+ # sort the current row by finding the element with the
+ # smallest absolute distance to each root in the
+ # previous row
+ available = range(len(prevrow))
+ for elem in row:
+ evect = elem-prevrow[available]
+ ind1 = abs(evect).argmin()
+ ind = available.pop(ind1)
+ sorted[n,ind] = elem
+ prevrow = sorted[n,:]
+ return sorted
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: Ryan K. <rk...@si...> - 2010-06-17 21:34:48
|
That sounds fine. I am totally open as far as licensing and what happens to the code. I don't know if there is an official statement in the code, but I would consider it BSD. I forgot to mention that I have my controls module on github: http://github.com/ryanGT/controls I don't know if that is helpful or not as far as collaboration goes. Ryan -- Ryan Krauss, Ph.D. Assistant Professor Mechanical Engineering Southern Illinois University Edwardsville On Thu, Jun 17, 2010 at 4:31 PM, Richard Murray <mu...@cd...> wrote: > Sounds good. I added a copy of your module into the subversion tree so that everyone can access a common copy. Let me know if you prefer it not live here (it won't be part of the official python-control distributions in its current location, but is accessible to people via a subversion checkout from SourceForce). > > -richard > > On 17 Jun 2010, at 14:18 , Ryan Krauss wrote: > >> That sounds great. It is good to know other people are working on >> this. I have had a few masters students ticker with it, but have been >> more or less been working on my own. >> >> I would be glad to be on both lists, but please use my list email: >> rya...@gm.... >> >> My TransferFunction class is basically a signal.ltisys with num and >> den as poly1d's (also from scipy). I mainly use the RootLocus and >> FreqResp methods along with the feedback function in my courses. >> >>> * You have implemented a couple of functions that I would love to pull into the python-control >>> library, including root locus plots and c2d support. Pulling those in would be a good way to see >>> how much synergy there is (based on how much rewrite is required). If this is OK by you, I >>> could give it a shot and let you know how it goes. >> >> The zoh c2d stuff took some work and I am happy with how it turned out >> (but I am not a digital control expert). So, this seems like a great >> first step. >> >>> * Going the other way, I've been trying to write the python-control code so that any object that >>> has the right attributes and/or member functions can be used. So, for example, the nyquist() >>> function should (eventually) work on any object for which there is a freqresp() member function >>> or a num and den attribute. >> >> I have tried to follow a similar approach where anything that has a >> num and den can be used in various parts of the code. >> >> For SISO systems, I think the code pretty much does what I would like >> it to. It could definitely use some well defined tests and better >> documentation and the scrutiny of more eyes on it, but I am more or >> less happy with it. >> >> The biggest need from my point of view is state-space stuff and >> especially MIMO. Actually, I don't even know if my code can elegantly >> handle SIMO. I need some SIMO tools for a project I am working on >> this summer at it is a little painful. >> >> I am certainly open to helping add new features. My rlocfind is >> slightly clumsy because I don't know who to click on a point on a plot >> and have IPython get the coordinates (without creating a wxPython gui >> just for the controls analysis or something). >> >> I will download your code in the next couple of days and give you my >> thoughts on how difficult it would be to combine the two code bases. >> >> -- >> Ryan Krauss, Ph.D. >> Assistant Professor >> Mechanical Engineering >> Southern Illinois University Edwardsville >> >> >> >> On Thu, Jun 17, 2010 at 3:43 PM, Richard Murray <mu...@cd...> wrote: >>> Hi Ryan. Good timing! I had just been trolling the web looking for others who are implementing control packages using python. Your name had already come up when I did a check a year or so back, but I was thinking that I would contact you to find out how things are going. And before I did, your message arrived! >>> >>> It would be great to find a way to combine efforts. I've got a couple of people who have expressed interest in adding to the functionality of the python-control package, so one possibility would be for us to take some of the functions that you have written and incorporate them into what we have so far. That would give us a sense of where the remaining gaps are and would allow us to look for some ways to continue to combine efforts. >>> >>> There are a couple of possible next steps: >>> >>> * You have implemented a couple of functions that I would love to pull into the python-control library, including root locus plots and c2d support. Pulling those in would be a good way to see how much synergy there is (based on how much rewrite is required). If this is OK by you, I could give it a shot and let you know how it goes. >>> >>> * Going the other way, I've been trying to write the python-control code so that any object that has the right attributes and/or member functions can be used. So, for example, the nyquist() function should (eventually) work on any object for which there is a freqresp() member function or a num and den attribute. It would be interesting to try this out and see whether we I can modify the python-control functions so that one of your TransferFunction objects could be passed to python-control functions. >>> >>> * Along these lines, if there are particular functions that you have been thinking about implementing next, please let us know so that we can either put them higher on the list or put them lower (if you will be working on them). We haven't yet done a prioritization, but my rough goal is to be able to run all of the examples in the controls book I wrote with Astrom (so that we can start using python in our classes as well), followed by some of the more advanced tools we use in a second quarter controls course (includes optimal control, trajectory generation, Kalman filter + extensions, etc). >>> >>> * Getting more active, if you would like to be signed up for either the python-control-announce list (low volume, summaries of when things are added) or the python-control-discuss (higher volume, used by people working on the project), just let me know. I'm cc'ing the latter list so that others know about the conversation. >>> >>> * Finally, if you actually have time to participate in the integration effort yourself, this would be wonderful. The "obvious" thing to do is to take the union of your TransferFunction class and the python-control TransferFunction class, then start working on generalizing various things to work for state space. For SISO this is trivial (signal.ltisys already keeps both num/den and ABCD representations for all objects), but gets tricker for things that might be MIMO. >>> >>> OK, that's probably enough for now -:). Let me know what you think about any and all of this. It would be great to collaborate! >>> >>> -richard >>> >>> On 17 Jun 2010, at 13:08 , Ryan Krauss wrote: >>> >>>> Richard, >>>> >>>> It looks like you are still actively developing a controls module for >>>> Python. I have used mine for teaching a first undergraduate controls >>>> course twice now and am relatively happy with it. It is poorly >>>> documented and doesn't really support state-space stuff yet, but it >>>> basically works for me. Let me know if you think it makes sense to >>>> combined our efforts somehow (I don't have a ton of time to spend on >>>> this, but would love to see other professors use Python for controls). >>>> >>>> Thanks, >>>> >>>> Ryan >>>> >>>> p.s. The version on my school webpage is slightly old. I need to fix >>>> that. Here are links to the temporary home of the latest versions >>>> from the course webpage: >>>> http://www.cs.siue.edu/~rkrauss/450/2010/controls_dist/controls-1.1.3.tar.gz >>>> http://www.cs.siue.edu/~rkrauss/450/2010/controls_dist/controls-1.1.3.win32.exe >>>> >>>> -- >>>> Ryan Krauss, Ph.D. >>>> Assistant Professor >>>> Mechanical Engineering >>>> Southern Illinois University Edwardsville >>> >>> > > |
|
From: Richard M. <mu...@cd...> - 2010-06-17 21:31:12
|
Sounds good. I added a copy of your module into the subversion tree so that everyone can access a common copy. Let me know if you prefer it not live here (it won't be part of the official python-control distributions in its current location, but is accessible to people via a subversion checkout from SourceForce). -richard On 17 Jun 2010, at 14:18 , Ryan Krauss wrote: > That sounds great. It is good to know other people are working on > this. I have had a few masters students ticker with it, but have been > more or less been working on my own. > > I would be glad to be on both lists, but please use my list email: > rya...@gm.... > > My TransferFunction class is basically a signal.ltisys with num and > den as poly1d's (also from scipy). I mainly use the RootLocus and > FreqResp methods along with the feedback function in my courses. > >> * You have implemented a couple of functions that I would love to pull into the python-control >> library, including root locus plots and c2d support. Pulling those in would be a good way to see >> how much synergy there is (based on how much rewrite is required). If this is OK by you, I >> could give it a shot and let you know how it goes. > > The zoh c2d stuff took some work and I am happy with how it turned out > (but I am not a digital control expert). So, this seems like a great > first step. > >> * Going the other way, I've been trying to write the python-control code so that any object that >> has the right attributes and/or member functions can be used. So, for example, the nyquist() >> function should (eventually) work on any object for which there is a freqresp() member function >> or a num and den attribute. > > I have tried to follow a similar approach where anything that has a > num and den can be used in various parts of the code. > > For SISO systems, I think the code pretty much does what I would like > it to. It could definitely use some well defined tests and better > documentation and the scrutiny of more eyes on it, but I am more or > less happy with it. > > The biggest need from my point of view is state-space stuff and > especially MIMO. Actually, I don't even know if my code can elegantly > handle SIMO. I need some SIMO tools for a project I am working on > this summer at it is a little painful. > > I am certainly open to helping add new features. My rlocfind is > slightly clumsy because I don't know who to click on a point on a plot > and have IPython get the coordinates (without creating a wxPython gui > just for the controls analysis or something). > > I will download your code in the next couple of days and give you my > thoughts on how difficult it would be to combine the two code bases. > > -- > Ryan Krauss, Ph.D. > Assistant Professor > Mechanical Engineering > Southern Illinois University Edwardsville > > > > On Thu, Jun 17, 2010 at 3:43 PM, Richard Murray <mu...@cd...> wrote: >> Hi Ryan. Good timing! I had just been trolling the web looking for others who are implementing control packages using python. Your name had already come up when I did a check a year or so back, but I was thinking that I would contact you to find out how things are going. And before I did, your message arrived! >> >> It would be great to find a way to combine efforts. I've got a couple of people who have expressed interest in adding to the functionality of the python-control package, so one possibility would be for us to take some of the functions that you have written and incorporate them into what we have so far. That would give us a sense of where the remaining gaps are and would allow us to look for some ways to continue to combine efforts. >> >> There are a couple of possible next steps: >> >> * You have implemented a couple of functions that I would love to pull into the python-control library, including root locus plots and c2d support. Pulling those in would be a good way to see how much synergy there is (based on how much rewrite is required). If this is OK by you, I could give it a shot and let you know how it goes. >> >> * Going the other way, I've been trying to write the python-control code so that any object that has the right attributes and/or member functions can be used. So, for example, the nyquist() function should (eventually) work on any object for which there is a freqresp() member function or a num and den attribute. It would be interesting to try this out and see whether we I can modify the python-control functions so that one of your TransferFunction objects could be passed to python-control functions. >> >> * Along these lines, if there are particular functions that you have been thinking about implementing next, please let us know so that we can either put them higher on the list or put them lower (if you will be working on them). We haven't yet done a prioritization, but my rough goal is to be able to run all of the examples in the controls book I wrote with Astrom (so that we can start using python in our classes as well), followed by some of the more advanced tools we use in a second quarter controls course (includes optimal control, trajectory generation, Kalman filter + extensions, etc). >> >> * Getting more active, if you would like to be signed up for either the python-control-announce list (low volume, summaries of when things are added) or the python-control-discuss (higher volume, used by people working on the project), just let me know. I'm cc'ing the latter list so that others know about the conversation. >> >> * Finally, if you actually have time to participate in the integration effort yourself, this would be wonderful. The "obvious" thing to do is to take the union of your TransferFunction class and the python-control TransferFunction class, then start working on generalizing various things to work for state space. For SISO this is trivial (signal.ltisys already keeps both num/den and ABCD representations for all objects), but gets tricker for things that might be MIMO. >> >> OK, that's probably enough for now -:). Let me know what you think about any and all of this. It would be great to collaborate! >> >> -richard >> >> On 17 Jun 2010, at 13:08 , Ryan Krauss wrote: >> >>> Richard, >>> >>> It looks like you are still actively developing a controls module for >>> Python. I have used mine for teaching a first undergraduate controls >>> course twice now and am relatively happy with it. It is poorly >>> documented and doesn't really support state-space stuff yet, but it >>> basically works for me. Let me know if you think it makes sense to >>> combined our efforts somehow (I don't have a ton of time to spend on >>> this, but would love to see other professors use Python for controls). >>> >>> Thanks, >>> >>> Ryan >>> >>> p.s. The version on my school webpage is slightly old. I need to fix >>> that. Here are links to the temporary home of the latest versions >>> from the course webpage: >>> http://www.cs.siue.edu/~rkrauss/450/2010/controls_dist/controls-1.1.3.tar.gz >>> http://www.cs.siue.edu/~rkrauss/450/2010/controls_dist/controls-1.1.3.win32.exe >>> >>> -- >>> Ryan Krauss, Ph.D. >>> Assistant Professor >>> Mechanical Engineering >>> Southern Illinois University Edwardsville >> >> |
|
From: <mur...@us...> - 2010-06-17 21:29:29
|
Revision: 26
http://python-control.svn.sourceforge.net/python-control/?rev=26&view=rev
Author: murrayrm
Date: 2010-06-17 21:29:23 +0000 (Thu, 17 Jun 2010)
Log Message:
-----------
added Ryan Krauss's controls.py module for easy access (and version control)
Added Paths:
-----------
trunk/external/
trunk/external/controls.py
Added: trunk/external/controls.py
===================================================================
--- trunk/external/controls.py (rev 0)
+++ trunk/external/controls.py 2010-06-17 21:29:23 UTC (rev 26)
@@ -0,0 +1,1386 @@
+# controls.py - Ryan Krauss's control module
+# $Id$
+
+"""This module is for analyzing linear, time-invariant dynamic systems
+and feedback control systems using the Laplace transform. The heart
+of the module is the TransferFunction class, which represents a
+transfer function as a ratio of numerator and denominator polynomials
+in s. TransferFunction is derived from scipy.signal.lti."""
+
+import glob, pdb
+from math import atan2, log10
+
+from scipy import *
+from scipy.linalg import inv as inverse
+from scipy.optimize import newton, fmin, fminbound
+from scipy.io import read_array, save, loadmat, write_array
+from scipy import signal
+
+from IPython.Debugger import Pdb
+
+import sys, os, copy, time
+
+from matplotlib.ticker import LogFormatterMathtext
+
+version = '1.1.0'
+
+class MyFormatter(LogFormatterMathtext):
+ def __call__(self, x, pos=None):
+ if pos==0: return '' # pos=0 is the first tick
+ else: return LogFormatterMathtext.__call__(self, x, pos)
+
+
+def shift(vectin, new):
+ N = len(vectin)-1
+ for n in range(N,0,-1):
+ vectin[n]=vectin[n-1]
+ vectin[0]=new
+ return vectin
+
+def myeq(p1,p2):
+ """Test the equality of the of two polynomials based on
+ coeffiecents."""
+ if hasattr(p1, 'coeffs') and hasattr(p2, 'coeffs'):
+ c1=p1.coeffs
+ c2=p2.coeffs
+ else:
+ return False
+ if len(c1)!=len(c2):
+ return False
+ else:
+ testvect=c1==c2
+ if hasattr(testvect,'all'):
+ return testvect.all()
+ else:
+ return testvect
+
+def build_fit_matrix(output_vect, input_vect, numorder, denorder):
+ """Build the [A] matrix used in least squares curve fitting
+ according to
+
+ output_vect = [A]c
+
+ as described in fit_discrete_response."""
+ A = zeros((len(output_vect),numorder+denorder+1))#the +1 accounts
+ #for the fact that both the numerator and the denominator
+ #have zero-order terms (which would give +2), but the
+ #zero order denominator term is actually not used in the fit
+ #(that is the output vector)
+ curin = input_vect
+ A[:,0] = curin
+ for n in range(1, numorder+1):
+ curin = r_[[0.0], curin[0:-1]]#prepend a 0 to curin and drop its
+ #last element
+ A[:,n] = curin
+ curout = -output_vect#this is the first output column, but it not
+ #actually used
+ firstden = numorder+1
+ for n in range(0, denorder):
+ curout = r_[[0.0], curout[0:-1]]
+ A[:,firstden+n] = curout
+ return A
+
+
+def fit_discrete_response(output_vect, input_vect, numorder, denorder):
+ """Find the coefficients of a digital transfer function that give
+ the best fit to output_vect in a least squares sense. output_vect
+ is the output of the system and input_vect is the input. The
+ input and output vectors are shifted backward in time a maximum of
+ numorder and denorder steps respectively. Each shifted vector
+ becomes a column in the matrix for the least squares curve fit of
+ the form
+
+ output_vect = [A]c
+
+ where [A] is the matrix whose columns are shifted versions of
+ input_vect and output_vect and c is composed of the numerator and
+ denominator coefficients of the transfer function. numorder and
+ denorder are the highest power of z in the numerator or
+ denominator respectively.
+
+ In essence, the approach is to find the coefficients that best fit
+ related the input_vect and output_vect according to the difference
+ equation
+
+ y(k) = b_0 x(k) + b_1 x(k-1) + b_2 x(k-2) + ... + b_m x(k-m)
+ - a_1 y(k-1) - a_2 y(k-2) - ... - a_n y(k-n)
+
+ where x = input_vect, y = output_vect, m = numorder, and
+ n = denorder. The unknown coefficient vector is then
+
+ c = [b_0, b_1, b_2, ... , b_m, a_1, a_2, ..., a_n]
+
+ Note that a_0 is forced to be 1.
+
+ The matrix [A] is then composed of [A] = [X(k) X(k-1) X(k-2)
+ ... Y(k-1) Y(k-2) ...] where X(k-2) represents the input_vect
+ shifted 2 elements and Y(k-2) represents the output_vect shifted
+ two elements."""
+ A = build_fit_matrix(output_vect, input_vect, numorder, denorder)
+ fitres = linalg.lstsq(A, output_vect)
+ x = fitres[0]
+ numz = x[0:numorder+1]
+ denz = x[numorder+1:]
+ denz = r_[[1.0],denz]
+ return numz, denz
+
+def prependzeros(num, den):
+ nd = len(den)
+ nn = len(num)
+ if nn < nd:
+ zvect = zeros(nd-nn)
+ numout = r_[zvect, num]
+ else:
+ numout = num
+ return numout, den
+
+def in_with_tol(elem, searchlist, rtol=1e-5, atol=1e-10):
+ """Determine whether or not elem+/-tol matches an element of
+ searchlist."""
+ for n, item in enumerate(searchlist):
+ if allclose(item, elem, rtol=rtol, atol=atol):
+ return n
+ return -1
+
+
+
+def PolyToLatex(polyin, var='s', fmt='%0.5g', eps=1e-12):
+ N = polyin.order
+ clist = polyin.coeffs
+ outstr = ''
+ for i, c in enumerate(clist):
+ curexp = N-i
+ curcoeff = fmt%c
+ if curexp > 0:
+ if curexp == 1:
+ curs = var
+ else:
+ curs = var+'^%i'%curexp
+ #Handle coeffs of +/- 1 in a special way:
+ if 1-eps < c < 1+eps:
+ curcoeff = ''
+ elif -1-eps < c < -1+eps:
+ curcoeff = '-'
+ else:
+ curs=''
+ curstr = curcoeff+curs
+ if c > 0 and outstr:
+ curcoeff = '+'+curcoeff
+ if abs(c) > eps:
+ outstr+=curcoeff+curs
+ return outstr
+
+
+def polyfactor(num, den, prepend=True, rtol=1e-5, atol=1e-10):
+ """Factor out any common roots from the polynomials represented by
+ the vectors num and den and return new coefficient vectors with
+ any common roots cancelled.
+
+ Because poly1d does not think in terms of z^-1, z^-2, etc. it may
+ be necessary to add zeros to the beginning of the numpoly coeffs
+ to represent multiplying through be z^-n where n is the order of
+ the denominator. If prependzeros is Trus, the numerator and
+ denominator coefficient vectors will have the same length."""
+ numpoly = poly1d(num)
+ denpoly = poly1d(den)
+ nroots = roots(numpoly).tolist()
+ droots = roots(denpoly).tolist()
+ n = 0
+ while n < len(nroots):
+ curn = nroots[n]
+ ind = in_with_tol(curn, droots, rtol=rtol, atol=atol)
+ if ind > -1:
+ nroots.pop(n)
+ droots.pop(ind)
+ #numpoly, rn = polydiv(numpoly, poly(curn))
+ #denpoly, rd = polydiv(denpoly, poly(curn))
+ else:
+ n += 1
+ numpoly = poly(nroots)
+ denpoly = poly(droots)
+ nvect = numpoly
+ dvect = denpoly
+ if prepend:
+ nout, dout = prependzeros(nvect, dvect)
+ else:
+ nout = nvect
+ dout = dvect
+ return nout, dout
+
+
+def polysubstitute(polyin, numsub, densub):
+ """Substitute one polynomial into another to support Tustin and
+ other c2d algorithms of a similar approach. The idea is to make
+ it easy to substitute
+
+ a z-1
+ s = - -----
+ T z+1
+
+ or other forms involving ratios of polynomials for s in a
+ polynomial of s such as the numerator or denominator of a transfer
+ function.
+
+ For the tustin example above, numsub=a*(z-1) and densub=T*(z+1),
+ where numsub and densub are scipy.poly1d instances.
+
+ Note that this approach seems to have substantial floating point
+ problems."""
+ mys = TransferFunction(numsub, densub)
+ out = 0.0
+ no = polyin.order
+ for n, coeff in enumerate(polyin.coeffs):
+ curterm = coeff*mys**(no-n)
+ out = out+curterm
+ return out
+
+
+def tustin_sub(polyin, T, a=2.0):
+ numsub = a*poly1d([1.0,-1.0])
+ densub = T*poly1d([1.0,1.0])
+ out = polysubstitute(polyin, numsub, densub)
+ out.myvar = 'z'
+ return out
+
+
+def create_swept_sine_input(maxt, dt, maxf, minf=0.0, deadtime=2.0):
+ t = arange(0, maxt, dt)
+ u = sweptsine(t, minf=minf, maxf=maxf)
+ if deadtime:
+ deadt = arange(0,deadtime, dt)
+ zv = zeros_like(deadt)
+ u = r_[zv, u, zv]
+ return u
+
+def create_swept_sine_t(maxt, dt, deadtime=2.0):
+ t = arange(0, maxt, dt)
+ if deadtime:
+ deadt = arange(0,deadtime, dt)
+ t = t+max(deadt)+dt
+ tpost = deadt+max(t)+dt
+ return r_[deadt, t, tpost]
+ else:
+ return t
+
+def ADC(vectin, bits=9, vmax=2.5, vmin=-2.5):
+ """Simulate the sampling portion of an analog-to-digital
+ conversion by outputing an integer number of counts associate with
+ each voltage in vectin."""
+ dv = (vmax-vmin)/2**bits
+ vect2 = clip(vectin, vmin, vmax)
+ counts = vect2/dv
+ return counts.astype(int)
+
+
+def CountsToFloat(counts, bits=9, vmax=2.5, vmin=-2.5):
+ """Convert the integer output of ADC to a floating point number by
+ mulitplying by dv."""
+ dv = (vmax-vmin)/2**bits
+ return dv*counts
+
+
+def epslist(listin, eps=1.0e-12):
+ """Make a copy of listin and then check each element of the copy
+ to see if its absolute value is greater than eps. Set to zero all
+ elements in the copied list whose absolute values are less than
+ eps. Return the copied list."""
+ listout = copy.deepcopy(listin)
+ for i in range(len(listout)):
+ if abs(listout[i])<eps:
+ listout[i] = 0.0
+ return listout
+
+
+def _PlotMatrixvsF(freqvect,matin,linetype='',linewidth=None, semilogx=True, allsolid=False, axis=None):
+ mykwargs={}
+ usepylab = False
+ if axis is None:
+ import pylab
+ axis = pylab.gca()
+ usepylab = True
+ if len(shape(matin))==1:
+ myargs=[freqvect,matin]
+ if linetype:
+ myargs.append(linetype)
+ else:
+ mykwargs.update(_getlinetype(axis))
+ if linewidth:
+ mykwargs['linewidth']=linewidth
+ if semilogx:
+ curline,=axis.semilogx(*myargs,**mykwargs)
+ else:
+ curline,=axis.plot(*myargs,**mykwargs)
+ mylines=[curline]
+# _inccount()
+ else:
+ mylines=[]
+ for q in range(shape(matin)[1]):
+ myargs=[freqvect,matin[:,q]]
+ if linetype:
+ myargs.append(linetype)
+ else:
+ mykwargs.update(_getlinetype(axis))
+ if linewidth:
+ mykwargs['linewidth']=linewidth
+ if semilogx:
+ curline,=axis.semilogx(*myargs,**mykwargs)
+ else:
+ curline,=axis.plot(*myargs,**mykwargs)
+ mylines.append(curline)
+# _inccount()
+ return mylines
+
+
+def _PlotMag(freqvect, bodein, linetype='-', linewidth=0, axis=None):
+ if callable(bodein.dBmag):
+ myvect=bodein.dBmag()
+ else:
+ myvect=bodein.dBmag
+ return _PlotMatrixvsF(freqvect, myvect, linetype=linetype, linewidth=linewidth, axis=axis)
+
+
+def _PlotPhase(freqvect, bodein, linetype='-', linewidth=0, axis=None):
+ return _PlotMatrixvsF(freqvect,bodein.phase,linetype=linetype,linewidth=linewidth, axis=axis)
+
+
+def _k_poles(TF,poleloc):
+ L = TF.num(poleloc)/TF.den(poleloc)
+ k = 1.0/abs(L)
+ poles = TF._RLFindRoots([k])
+ poles = TF._RLSortRoots(poles)
+ return k,poles
+
+def _checkpoles(poleloc,pnew):
+ evect = abs(poleloc-array(pnew))
+ ind = evect.argmin()
+ pout = pnew[ind]
+ return pout
+
+
+def _realizable(num, den):
+ realizable = False
+ if not isscalar(den):
+ if isscalar(num):
+ realizable = True
+ elif len(den) >= len(num):
+ realizable = True
+ return realizable
+
+
+def shape_u(uvect, slope):
+ u_shaped = zeros_like(uvect)
+ u_shaped[0] = uvect[0]
+
+ N = len(uvect)
+
+ for n in range(1, N):
+ diff = uvect[n] - u_shaped[n-1]
+ if diff > slope:
+ u_shaped[n] = u_shaped[n-1] + slope
+ elif diff < -1*slope:
+ u_shaped[n] = u_shaped[n-1] - slope
+ else:
+ u_shaped[n] = uvect[n]
+ return u_shaped
+
+
+class TransferFunction(signal.lti):
+ def __setattr__(self, attr, val):
+ realizable = False
+ if hasattr(self, 'den') and hasattr(self, 'num'):
+ realizable = _realizable(self.num, self.den)
+ if realizable:
+ signal.lti.__setattr__(self, attr, val)
+ else:
+ self.__dict__[attr] = val
+
+
+ def __init__(self, num, den, dt=0.01, maxt=5.0, myvar='s'):
+ """num and den are either scalar constants or lists that are
+ passed to scipy.poly1d to create a list of coefficients."""
+ #print('in TransferFunction.__init__, dt=%s' % dt)
+ if _realizable(num, den):
+ signal.lti.__init__(self, num, den)
+ self.num = poly1d(num)
+ self.den = poly1d(den)
+ self.dt = dt
+ self.myvar = myvar
+ self.maxt = maxt
+
+
+ def __repr__(self, labelstr='controls.TransferFunction'):
+ nstr=str(self.num)#.strip()
+ dstr=str(self.den)#.strip()
+ nstr=nstr.replace('x',self.myvar)
+ dstr=dstr.replace('x',self.myvar)
+ n=len(dstr)
+ m=len(nstr)
+ shift=(n-m)/2*' '
+ nstr=nstr.replace('\n','\n'+shift)
+ tempstr=labelstr+'\n'+shift+nstr+'\n'+'-'*n+'\n '+dstr
+ return tempstr
+
+
+ def __call__(self,s,optargs=()):
+ return self.num(s)/self.den(s)
+
+
+ def __add__(self,other):
+ if hasattr(other,'num') and hasattr(other,'den'):
+ if len(self.den.coeffs)==len(other.den.coeffs) and \
+ (self.den.coeffs==other.den.coeffs).all():
+ return TransferFunction(self.num+other.num,self.den)
+ else:
+ return TransferFunction(self.num*other.den+other.num*self.den,self.den*other.den)
+ elif isinstance(other, int) or isinstance(other, float):
+ return TransferFunction(other*self.den+self.num,self.den)
+ else:
+ raise ValueError, 'do not know how to add TransferFunction and '+str(other) +' which is of type '+str(type(other))
+
+ def __radd__(self,other):
+ return self.__add__(other)
+
+
+ def __mul__(self,other):
+ if isinstance(other, Digital_P_Control):
+ return self.__class__(other.kp*self.num, self.den)
+ elif hasattr(other,'num') and hasattr(other,'den'):
+ if myeq(self.num,other.den) and myeq(self.den,other.num):
+ return 1
+ elif myeq(self.num,other.den):
+ return self.__class__(other.num,self.den)
+ elif myeq(self.den,other.num):
+ return self.__class__(self.num,other.den)
+ else:
+ gain = self.gain*other.gain
+ new_num, new_den = polyfactor(self.num*other.num, \
+ self.den*other.den)
+ newtf = self.__class__(new_num*gain, new_den)
+ return newtf
+ elif isinstance(other, int) or isinstance(other, float):
+ return self.__class__(other*self.num,self.den)
+
+
+ def __pow__(self, expon):
+ """Basically, go self*self*self as many times as necessary. I
+ haven't thought about negative exponents. I don't think this
+ would be hard, you would just need to keep dividing by self
+ until you got the right answer."""
+ assert expon >= 0, 'TransferFunction.__pow__ does not yet support negative exponents.'
+ out = 1.0
+ for n in range(expon):
+ out *= self
+ return out
+
+
+ def __rmul__(self,other):
+ return self.__mul__(other)
+
+
+ def __div__(self,other):
+ if hasattr(other,'num') and hasattr(other,'den'):
+ if myeq(self.den,other.den):
+ return TransferFunction(self.num,other.num)
+ else:
+ return TransferFunction(self.num*other.den,self.den*other.num)
+ elif isinstance(other, int) or isinstance(other, float):
+ return TransferFunction(self.num,other*self.den)
+
+
+ def __rdiv__(self, other):
+ print('calling TransferFunction.__rdiv__')
+ return self.__div__(other)
+
+
+ def __truediv__(self,other):
+ return self.__div__(other)
+
+
+ def _get_set_dt(self, dt=None):
+ if dt is not None:
+ self.dt = float(dt)
+ return self.dt
+
+
+ def ToLatex(self, eps=1e-12, fmt='%0.5g', ds=True):
+ mynum = self.num
+ myden = self.den
+ npart = PolyToLatex(mynum)
+ dpart = PolyToLatex(myden)
+ outstr = '\\frac{'+npart+'}{'+dpart+'}'
+ if ds:
+ outstr = '\\displaystyle '+outstr
+ return outstr
+
+
+ def RootLocus(self, kvect, fig=None, fignum=1, \
+ clear=True, xlim=None, ylim=None, plotstr='-'):
+ """Calculate the root locus by finding the roots of 1+k*TF(s)
+ where TF is self.num(s)/self.den(s) and each k is an element
+ of kvect."""
+ if fig is None:
+ import pylab
+ fig = pylab.figure(fignum)
+ if clear:
+ fig.clf()
+ ax = fig.add_subplot(111)
+ mymat = self._RLFindRoots(kvect)
+ mymat = self._RLSortRoots(mymat)
+ #plot open loop poles
+ poles = array(self.den.r)
+ ax.plot(real(poles), imag(poles), 'x')
+ #plot open loop zeros
+ zeros = array(self.num.r)
+ if zeros.any():
+ ax.plot(real(zeros), imag(zeros), 'o')
+ for col in mymat.T:
+ ax.plot(real(col), imag(col), plotstr)
+ if xlim:
+ ax.set_xlim(xlim)
+ if ylim:
+ ax.set_ylim(ylim)
+ ax.set_xlabel('Real')
+ ax.set_ylabel('Imaginary')
+ return mymat
+
+
+ def _RLFindRoots(self, kvect):
+ """Find the roots for the root locus."""
+ roots = []
+ for k in kvect:
+ curpoly = self.den+k*self.num
+ curroots = curpoly.r
+ curroots.sort()
+ roots.append(curroots)
+ mymat = row_stack(roots)
+ return mymat
+
+
+ def _RLSortRoots(self, mymat):
+ """Sort the roots from self._RLFindRoots, so that the root
+ locus doesn't show weird pseudo-branches as roots jump from
+ one branch to another."""
+ sorted = zeros_like(mymat)
+ for n, row in enumerate(mymat):
+ if n==0:
+ sorted[n,:] = row
+ else:
+ #sort the current row by finding the element with the
+ #smallest absolute distance to each root in the
+ #previous row
+ available = range(len(prevrow))
+ for elem in row:
+ evect = elem-prevrow[available]
+ ind1 = abs(evect).argmin()
+ ind = available.pop(ind1)
+ sorted[n,ind] = elem
+ prevrow = sorted[n,:]
+ return sorted
+
+
+ def opt(self, kguess):
+ pnew = self._RLFindRoots(kguess)
+ pnew = self._RLSortRoots(pnew)[0]
+ if len(pnew)>1:
+ pnew = _checkpoles(self.poleloc,pnew)
+ e = abs(pnew-self.poleloc)**2
+ return sum(e)
+
+
+ def rlocfind(self, poleloc):
+ self.poleloc = poleloc
+ kinit,pinit = _k_poles(self,poleloc)
+ k = optimize.fmin(self.opt,[kinit])[0]
+ poles = self._RLFindRoots([k])
+ poles = self._RLSortRoots(poles)
+ return k, poles
+
+
+ def PlotTimeResp(self, u, t, fig, clear=True, label='model', mysub=111):
+ ax = fig.add_subplot(mysub)
+ if clear:
+ ax.cla()
+ try:
+ y = self.lsim(u, t)
+ except:
+ y = self.lsim2(u, t)
+ ax.plot(t, y, label=label)
+ return ax
+
+
+## def BodePlot(self, f, fig, clear=False):
+## mtf = self.FreqResp(
+## ax1 = fig.axes[0]
+## ax1.semilogx(modelf,20*log10(abs(mtf)))
+## mphase = angle(mtf, deg=1)
+## ax2 = fig.axes[1]
+## ax2.semilogx(modelf, mphase)
+
+
+ def SimpleFactor(self):
+ mynum=self.num
+ myden=self.den
+ dsf=myden[myden.order]
+ nsf=mynum[mynum.order]
+ sden=myden/dsf
+ snum=mynum/nsf
+ poles=sden.r
+ residues=zeros(shape(sden.r),'D')
+ factors=[]
+ for x,p in enumerate(poles):
+ polearray=poles.copy()
+ polelist=polearray.tolist()
+ mypole=polelist.pop(x)
+ tempden=1.0
+ for cp in polelist:
+ tempden=tempden*(poly1d([1,-cp]))
+ tempTF=TransferFunction(snum,tempden)
+ curres=tempTF(mypole)
+ residues[x]=curres
+ curTF=TransferFunction(curres,poly1d([1,-mypole]))
+ factors.append(curTF)
+ return factors,nsf,dsf
+
+ def factor_constant(self, const):
+ """Divide numerator and denominator coefficients by const"""
+ self.num = self.num/const
+ self.den = self.den/const
+
+ def lsim(self, u, t, interp=0, returnall=False, X0=None):
+ """Find the response of the TransferFunction to the input u
+ with time vector t. Uses signal.lsim.
+
+ return y the response of the system."""
+ if returnall:#most users will just want the system output y,
+ #but some will need the (t, y, x) tuple that
+ #signal.lsim returns
+ return signal.lsim(self, u, t, interp=interp, X0=X0)
+ else:
+ return signal.lsim(self, u, t, interp=interp, X0=X0)[1]
+
+ def lsim2(self, u, t, returnall=False, X0=None):
+ #tempsys=signal.lti(self.num,self.den)
+ if returnall:
+ return signal.lsim2(self, u, t, X0=X0)
+ else:
+ return signal.lsim2(self, u, t, X0=X0)[1]
+
+
+ def residue(self, tol=1e-3, verbose=0):
+ """from scipy.signal.residue:
+
+ Compute residues/partial-fraction expansion of b(s) / a(s).
+
+ If M = len(b) and N = len(a)
+
+ b(s) b[0] s**(M-1) + b[1] s**(M-2) + ... + b[M-1]
+ H(s) = ------ = ----------------------------------------------
+ a(s) a[0] s**(N-1) + a[1] s**(N-2) + ... + a[N-1]
+
+ r[0] r[1] r[-1]
+ = -------- + -------- + ... + --------- + k(s)
+ (s-p[0]) (s-p[1]) (s-p[-1])
+
+ If there are any repeated roots (closer than tol), then the
+ partial fraction expansion has terms like
+
+ r[i] r[i+1] r[i+n-1]
+ -------- + ----------- + ... + -----------
+ (s-p[i]) (s-p[i])**2 (s-p[i])**n
+
+ returns r, p, k
+ """
+ r,p,k = signal.residue(self.num, self.den, tol=tol)
+ if verbose>0:
+ print('r='+str(r))
+ print('')
+ print('p='+str(p))
+ print('')
+ print('k='+str(k))
+
+ return r, p, k
+
+
+ def PartFrac(self, eps=1.0e-12):
+ """Compute the partial fraction expansion based on the residue
+ command. In the final polynomials, coefficients whose
+ absolute values are less than eps are set to zero."""
+ r,p,k = self.residue()
+
+ rlist = r.tolist()
+ plist = p.tolist()
+
+ N = len(rlist)
+
+ tflist = []
+ eps = 1e-12
+
+ while N > 0:
+ curr = rlist.pop(0)
+ curp = plist.pop(0)
+ if abs(curp.imag) < eps:
+ #This is a purely real pole. The portion of the partial
+ #fraction expansion corresponding to this pole is curr/(s-curp)
+ curtf = TransferFunction(curr,[1,-curp])
+ else:
+ #this is a complex pole and we need to find its conjugate and
+ #handle them together
+ cind = plist.index(curp.conjugate())
+ rconj = rlist.pop(cind)
+ pconj = plist.pop(cind)
+ p1 = poly1d([1,-curp])
+ p2 = poly1d([1,-pconj])
+ #num = curr*p2+rconj*p1
+ Nr = curr.real
+ Ni = curr.imag
+ Pr = curp.real
+ Pi = curp.imag
+ numlist = [2.0*Nr,-2.0*(Nr*Pr+Ni*Pi)]
+ numlist = epslist(numlist, eps)
+ num = poly1d(numlist)
+ denlist = [1, -2.0*Pr,Pr**2+Pi**2]
+ denlist = epslist(denlist, eps)
+ den = poly1d(denlist)
+ curtf = TransferFunction(num,den)
+ tflist.append(curtf)
+ N = len(rlist)
+ return tflist
+
+
+ def FreqResp(self, f, fignum=1, fig=None, clear=True, \
+ grid=True, legend=None, legloc=1, legsub=1, **kwargs):
+ """Compute the frequency response of the transfer function
+ using the frequency vector f, returning a complex vector.
+
+ The frequency response (Bode plot) will be plotted on
+ figure(fignum) unless fignum=None.
+
+ legend should be a list of legend entries if a legend is
+ desired. If legend is not None, the legend will be placed on
+ the top half of the plot (magnitude portion) if legsub=1, or
+ on the bottom half with legsub=2. legloc follows the same
+ rules as the pylab legend command (1 is top right and goes
+ counter-clockwise from there.)"""
+ testvect=real(f)==0
+ if testvect.all():
+ s=f#then you really sent me s and not f
+ else:
+ s=2.0j*pi*f
+ self.comp = self.num(s)/self.den(s)
+ self.dBmag = 20*log10(abs(self.comp))
+ self.phase = angle(self.comp,1)
+
+ if fig is None:
+ if fignum is not None:
+ import pylab
+ fig = pylab.figure(fignum)
+
+ if fig is not None:
+ if clear:
+ fig.clf()
+
+ if fig is not None:
+ myargs=['linetype','colors','linewidth']
+ subkwargs={}
+ for key in myargs:
+ if kwargs.has_key(key):
+ subkwargs[key]=kwargs[key]
+ if clear:
+ fig.clf()
+ ax1 = fig.add_subplot(2,1,1)
+ #if clear:
+ # ax1.cla()
+ myind=ax1._get_lines.count
+ mylines=_PlotMag(f, self, axis=ax1, **subkwargs)
+ ax1.set_ylabel('Mag. Ratio (dB)')
+ ax1.xaxis.set_major_formatter(MyFormatter())
+ if grid:
+ ax1.grid(1)
+ if legend is not None and legsub==1:
+ ax1.legend(legend, legloc)
+ ax2 = fig.add_subplot(2,1,2, sharex=ax1)
+ #if clear:
+ # ax2.cla()
+ mylines=_PlotPhase(f, self, axis=ax2, **subkwargs)
+ ax2.set_ylabel('Phase (deg.)')
+ ax2.set_xlabel('Freq. (Hz)')
+ ax2.xaxis.set_major_formatter(MyFormatter())
+ if grid:
+ ax2.grid(1)
+ if legend is not None and legsub==2:
+ ax2.legend(legend, legloc)
+ return self.comp
+
+
+ def CrossoverFreq(self, f):
+ if not hasattr(self, 'dBmag'):
+ self.FreqResp(f, fignum=None)
+ t1 = squeeze(self.dBmag > 0.0)
+ t2 = r_[t1[1:],t1[0]]
+ t3 = (t1 & -t2)
+ myinds = where(t3)[0]
+ if not myinds.any():
+ return None, []
+ maxind = max(myinds)
+ return f[maxind], maxind
+
+
+ def PhaseMargin(self,f):
+ fc,ind=self.CrossoverFreq(f)
+ if not fc:
+ return 180.0
+ return 180.0+squeeze(self.phase[ind])
+
+
+ def create_tvect(self, dt=None, maxt=None):
+ if dt is None:
+ dt = self.dt
+ else:
+ self.dt = dt
+ assert dt is not None, "You must either pass in a dt or call create_tvect on an instance with a self.dt already defined."
+ if maxt is None:
+ if hasattr(self,'maxt'):
+ maxt = self.maxt
+ else:
+ maxt = 100*dt
+ else:
+ self.maxt = maxt
+ tvect = arange(0,maxt+dt/2.0,dt)
+ self.t = tvect
+ return tvect
+
+
+ def create_impulse(self, dt=None, maxt=None, imp_time=0.5):
+ """Create the input impulse vector to be used in least squares
+ curve fitting of the c2d function."""
+ if dt is None:
+ dt = self.dt
+ indon = int(imp_time/dt)
+ tvect = self.create_tvect(dt=dt, maxt=maxt)
+ imp = zeros_like(tvect)
+ imp[indon] = 1.0
+ return imp
+
+
+ def create_step_input(self, dt=None, maxt=None, indon=5):
+ """Create the input impulse vector to be used in least squares
+ curve fitting of the c2d function."""
+ tvect = self.create_tvect(dt=dt, maxt=maxt)
+ mystep = zeros_like(tvect)
+ mystep[indon:] = 1.0
+ return mystep
+
+
+ def step_response(self, t=None, dt=None, maxt=None, \
+ step_time=None, fignum=1, clear=True, \
+ plotu=False, amp=1.0, interp=0, fig=None, \
+ fmts=['-','-'], legloc=0, returnall=0, \
+ legend=None, **kwargs):
+ """Find the response of the system to a step input. If t is
+ not given, then the time vector will go from 0 to maxt in
+ steps of dt i.e. t=arange(0,maxt,dt). If dt and maxt are not
+ given, the parameters from the TransferFunction instance will
+ be used.
+
+ step_time is the time when the step input turns on. If not
+ given, it will default to 0.
+
+ If clear is True, the figure will be cleared first.
+ clear=False could be used to overlay the step responses of
+ multiple TransferFunction's.
+
+ plotu=True means that the step input will also be shown on the
+ graph.
+
+ amp is the amplitude of the step input.
+
+ return y unless returnall is set then return y, t, u
+
+ where y is the response of the transfer function, t is the
+ time vector, and u is the step input vector."""
+ if t is not None:
+ tvect = t
+ else:
+ tvect = self.create_tvect(dt=dt, maxt=maxt)
+ u = zeros_like(tvect)
+ if dt is None:
+ dt = self.dt
+ if step_time is None:
+ step_time = 0.0
+ #step_time = 0.1*tvect.max()
+ if kwargs.has_key('indon'):
+ indon = kwargs['indon']
+ else:
+ indon = int(step_time/dt)
+ u[indon:] = amp
+ try:
+ ystep = self.lsim(u, tvect, interp=interp)#[1]#the outputs of lsim are (t, y,x)
+ except:
+ ystep = self.lsim2(u, tvect)#[1]
+
+ if fig is None:
+ if fignum is not None:
+ import pylab
+ fig = pylab.figure(fignum)
+
+ if fig is not None:
+ if clear:
+ fig.clf()
+ ax = fig.add_subplot(111)
+ if plotu:
+ leglist =['Input','Output']
+ ax.plot(tvect, u, fmts[0], linestyle='steps', **kwargs)#assume step input wants 'steps' linestyle
+ ofmt = fmts[1]
+ else:
+ ofmt = fmts[0]
+ ax.plot(tvect, ystep, ofmt, **kwargs)
+ ax.set_ylabel('Step Response')
+ ax.set_xlabel('Time (sec)')
+ if legend is not None:
+ ax.legend(legend, loc=legloc)
+ elif plotu:
+ ax.legend(leglist, loc=legloc)
+ #return ystep, ax
+ #else:
+ #return ystep
+ if returnall:
+ return ystep, tvect, u
+ else:
+ return ystep
+
+
+
+ def impulse_response(self, dt=None, maxt=None, fignum=1, \
+ clear=True, amp=1.0, fig=None, \
+ fmt='-', **kwargs):
+ """Find the impulse response of the system using
+ scipy.signal.impulse.
+
+ The time vector will go from 0 to maxt in steps of dt
+ i.e. t=arange(0,maxt,dt). If dt and maxt are not given, the
+ parameters from the TransferFunction instance will be used.
+
+ If clear is True, the figure will be cleared first.
+ clear=False could be used to overlay the impulse responses of
+ multiple TransferFunction's.
+
+ amp is the amplitude of the impulse input.
+
+ return y, t
+
+ where y is the impulse response of the transfer function and t
+ is the time vector."""
+
+ tvect = self.create_tvect(dt=dt, maxt=maxt)
+ temptf = amp*self
+ tout, yout = temptf.impulse(T=tvect)
+
+ if fig is None:
+ if fignum is not None:
+ import pylab
+ fig = pylab.figure(fignum)
+
+ if fig is not None:
+ if clear:
+ fig.clf()
+ ax = fig.add_subplot(111)
+ ax.plot(tvect, yout, fmt, **kwargs)
+ ax.set_ylabel('Impulse Response')
+ ax.set_xlabel('Time (sec)')
+
+ return yout, tout
+
+
+ def swept_sine_response(self, maxf, minf=0.0, dt=None, maxt=None, deadtime=2.0, interp=0):
+ u = create_swept_sine_input(maxt, dt, maxf, minf=minf, deadtime=deadtime)
+ t = create_swept_sine_t(maxt, dt, deadtime=deadtime)
+ ysweep = self.lsim(u, t, interp=interp)
+ return t, u, ysweep
+
+
+ def _c2d_sub(self, numsub, densub, scale):
+ """This method performs substitutions for continuous to
+ digital conversions using the form:
+
+ numsub
+ s = scale* --------
+ densub
+
+ where scale is a floating point number and numsub and densub
+ are poly1d instances.
+
+ For example, scale = 2.0/T, numsub = poly1d([1,-1]), and
+ densub = poly1d([1,1]) for a Tustin c2d transformation."""
+ m = self.num.order
+ n = self.den.order
+ mynum = 0.0
+ for p, coeff in enumerate(self.num.coeffs):
+ mynum += poly1d(coeff*(scale**(m-p))*((numsub**(m-p))*(densub**(n-(m-p)))))
+ myden = 0.0
+ for p, coeff in enumerate(self.den.coeffs):
+ myden += poly1d(coeff*(scale**(n-p))*((numsub**(n-p))*(densub**(n-(n-p)))))
+ return mynum.coeffs, myden.coeffs
+
+
+ def c2d_tustin(self, dt=None, a=2.0):
+ """Convert a continuous time transfer function into a digital
+ one by substituting
+
+ a z-1
+ s = - -----
+ T z+1
+
+ into the compensator, where a is typically 2.0"""
+ #print('in TransferFunction.c2d_tustin, dt=%s' % dt)
+ dt = self._get_set_dt(dt)
+ #print('in TransferFunction.c2d_tustin after _get_set_dt, dt=%s' % dt)
+ scale = a/dt
+ numsub = poly1d([1.0,-1.0])
+ densub = poly1d([1.0,1.0])
+ mynum, myden = self._c2d_sub(numsub, densub, scale)
+ mynum = mynum/myden[0]
+ myden = myden/myden[0]
+ return mynum, myden
+
+
+
+ def c2d(self, dt=None, maxt=None, method='zoh', step_time=0.5, a=2.0):
+ """Find a numeric approximation of the discrete transfer
+ function of the system.
+
+ The general approach is to find the response of the system
+ using lsim and fit a discrete transfer function to that
+ response as a least squares problem.
+
+ dt is the time between discrete time intervals (i.e. the
+ sample time).
+
+ maxt is the length of time for which to calculate the system
+ respnose. An attempt is made to guess an appropriate stopping
+ time if maxt is None. For now, this defaults to 100*dt,
+ assuming that dt is appropriate for the system poles.
+
+ method is a string describing the c2d conversion algorithm.
+ method = 'zoh refers to a zero-order hold for a sampled-data
+ system and follows the approach outlined by Dorsey in section
+ 14.17 of
+ "Continuous and Discrete Control Systems" summarized on page
+ 472 of the 2002 edition.
+
+ Other supported options for method include 'tustin'
+
+ indon gives the index of when the step input should switch on
+ for zoh or when the impulse should happen otherwise. There
+ should probably be enough zero entries before the input occurs
+ to accomidate the order of the discrete transfer function.
+
+ a is used only if method = 'tustin' and it is substituted in the form
+
+ a z-1
+ s = - -----
+ T z+1
+
+ a is almost always equal to 2.
+ """
+ if method.lower() == 'zoh':
+ ystep = self.step_response(dt=dt, maxt=maxt, step_time=step_time)[0]
+ myimp = self.create_impulse(dt=dt, maxt=maxt, imp_time=step_time)
+ #Pdb().set_trace()
+ print('You called c2d with "zoh". This is most likely bad.')
+ nz, dz = fit_discrete_response(ystep, myimp, self.den.order, self.den.order+1)#we want the numerator order to be one less than the denominator order - the denominator order +1 is the order of the denominator during a step response
+ #multiply by (1-z^-1)
+ nz2 = r_[nz, [0.0]]
+ nzs = r_[[0.0],nz]
+ nz3 = nz2 - nzs
+ nzout, dzout = polyfactor(nz3, dz)
+ return nzout, dzout
+ #return nz3, dz
+ elif method.lower() == 'tustin':
+ #The basic approach for tustin is to create a transfer
+ #function that represents s mapped into z and then
+ #substitute this s(z)=a/T*(z-1)/(z+1) into the continuous
+ #transfer function
+ return self.c2d_tustin(dt=dt, a=a)
+ else:
+ raise ValueError, 'c2d method not understood:'+str(method)
+
+
+
+ def DigitalSim(self, u, method='zoh', bits=9, vmin=-2.5, vmax=2.5, dt=None, maxt=None, digitize=True):
+ """Simulate the digital reponse of the transfer to input u. u
+ is assumed to be an input signal that has been sampled with
+ frequency 1/dt. u is further assumed to be a floating point
+ number with precision much higher than bits. u will be
+ digitized over the range [min, max], which is broken up into
+ 2**bits number of bins.
+
+ The A and B vectors from c2d conversion will be found using
+ method, dt, and maxt. Note that maxt is only used for
+ method='zoh'.
+
+ Once A and B have been found, the digital reponse of the
+ system to the digitized input u will be found."""
+ B, A = self.c2d(dt=dt, maxt=maxt, method=method)
+ assert A[0]==1.0, "A[0]!=1 in c2d result, A="+str(A)
+ uvect = zeros(len(B), dtype='d')
+ yvect = zeros(len(A)-1, dtype='d')
+ if digitize:
+ udig = ADC(u, bits, vmax=vmax, vmin=vmin)
+ dv = (vmax-vmin)/(2**bits-1)
+ else:
+ udig = u
+ dv = 1.0
+ Ydig = zeros(len(u), dtype='d')
+ for n, u0 in enumerate(udig):
+ uvect = shift(uvect, u0)
+ curY = dot(uvect,B)
+ negpart = dot(yvect,A[1:])
+ curY -= negpart
+ if digitize:
+ curY = int(curY)
+ Ydig[n] = curY
+ yvect = shift(yvect, curY)
+ return Ydig*dv
+
+
+
+class Input(TransferFunction):
+ def __repr__(self):
+ return TransferFunction.__repr__(self, labelstr='controls.Input')
+
+
+class Compensator(TransferFunction):
+ def __init__(self, num, den, *args, **kwargs):
+ #print('in Compensator.__init__')
+ #Pdb().set_trace()
+ TransferFunction.__init__(self, num, den, *args, **kwargs)
+
+
+ def c2d(self, dt=None, a=2.0):
+ """Compensators should use Tustin for c2d conversion. This
+ method is just and alias for TransferFunction.c2d_tustin"""
+ #print('in Compensators.c2d, dt=%s' % dt)
+ #Pdb().set_trace()
+ return TransferFunction.c2d_tustin(self, dt=dt, a=a)
+
+ def __repr__(self):
+ return TransferFunction.__repr__(self, labelstr='controls.Compensator')
+
+
+
+class Digital_Compensator(object):
+ def __init__(self, num, den, input_vect=None, output_vect=None):
+ self.num = num
+ self.den = den
+ self.input = input_vect
+ self.output = output_vect
+ self.Nnum = len(self.num)
+ self.Nden = len(self.den)
+
+
+ def calc_out(self, i):
+ out = 0.0
+ for n, bn in enumerate(self.num):
+ out += self.input[i-n]*bn
+
+ for n in range(1, self.Nden):
+ out -= self.output[i-n]*self.den[n]
+ out = out/self.den[0]
+ return out
+
+
+class Digital_PI(object):
+ def __init__(self, kp, ki, input_vect=None, output_vect=None):
+ self.kp = kp
+ self.ki = ki
+ self.input = input_vect
+ self.output = output_vect
+ self.esum = 0.0
+
+
+ def prep(self):
+ self.esum = zeros_like(self.input)
+
+
+ def calc_out(self, i):
+ self.esum[i] = self.esum[i-1]+self.input[i]
+ out = self.input[i]*self.kp+self.esum[i]*self.ki
+ return out
+
+
+class Digital_P_Control(Digital_Compensator):
+ def __init__(self, kp, input_vect=None, output_vect=None):
+ self.kp = kp
+ self.input = input_vect
+ self.output = output_vect
+ self.num = poly1d([kp])
+ self.den = poly1d([1])
+ self.gain = 1
+
+ def calc_out(self, i):
+ self.output[i] = self.kp*self.input[i]
+ return self.output[i]
+
+
+def dig_comp_from_c_comp(c_comp, dt):
+ """Convert a continuous compensator into a digital one using Tustin
+ and sampling time dt."""
+ b, a = c_comp.c2d_tustin(dt=dt)
+ return Digital_Compensator(b, a)
+
+
+class FirstOrderCompensator(Compensator):
+ def __init__(self, K, z, p, dt=0.004):
+ """Create a first order compensator whose transfer function is
+
+ K*(s+z)
+ D(s) = -----------
+ (s+p) """
+ Compensator.__init__(self, K*poly1d([1,z]), [1,p])
+
+
+ def __repr__(self):
+ return TransferFunction.__repr__(self, labelstr='controls.FirstOrderCompensator')
+
+
+ def ToPSoC(self, dt=0.004):
+ b, a = self.c2d(dt=dt)
+ outstr = 'v = %f*e%+f*ep%+f*vp;'%(b[0],b[1],-a[1])
+ print('PSoC str:')
+ print(outstr)
+ return outstr
+
+
+def sat(vin, vmax=2.0):
+ if vin > vmax:
+ return vmax
+ elif vin < -1*vmax:
+ return -1*vmax
+ else:
+ return vin
+
+
+
+class Closed_Loop_System_with_Sat(object):
+ def __init__(self, plant_tf, Kp, sat):
+ self.plant_tf = plant_tf
+ self.Kp = Kp
+ self.sat = sat
+
+
+ def lsim(self, u, t, X0=None, include_sat=True, \
+ returnall=0, lsim2=0, verbosity=0):
+ dt = t[1]-t[0]
+ if X0 is None:
+ X0 = zeros((2,len(self.plant_tf.den.coeffs)-1))
+ N = len(t)
+ y = zeros(N)
+ v = zeros(N)
+ x_n = X0
+ for n in range(1,N):
+ t_n = t[n]
+ if verbosity > 0:
+ print('t_n='+str(t_n))
+ e = u[n]-y[n-1]
+ v_n = self.Kp*e
+ if include_sat:
+ v_n = sat(v_n, vmax=self.sat)
+ #simulate for one dt using ZOH
+ if lsim2:
+ t_nn, y_n, x_n = self.plant_tf.lsim2([v_n,v_n], [t_n, t_n+dt], X0=x_n[-1], returnall=1)
+ else:
+ t_nn, y_n, x_n = self.plant_tf.lsim([v_n,v_n], [t_n, t_n+dt], X0=x_n[-1], returnall=1)
+
+ y[n] = y_n[-1]
+ v[n] = v_n
+ self.y = y
+ self.v = v
+ self.u = u
+ if returnall:
+ return y, v
+ else:
+ return y
+
+
+
+
+
+def step_input():
+ return Input(1,[1,0])
+
+
+def feedback(olsys,H=1):
+ """Calculate the closed-loop transfer function
+
+ olsys
+ cltf = --------------
+ 1+H*olsys
+
+ where olsys is the transfer function of the open loop
+ system (Gc*Gp) and H is the transfer function in the feedback
+ loop (H=1 for unity feedback)."""
+ clsys=olsys/(1.0+H*olsys)
+ return clsys
+
+
+
+def Usweep(ti,maxt,minf=0.0,maxf=10.0):
+ """Return the current value (scalar) of a swept sine signal - must be used
+ with list comprehension to generate a vector.
+
+ ti - current time (scalar)
+ minf - lowest frequency in the sweep
+ maxf - highest frequency in the sweep
+ maxt - T or the highest value in the time vector"""
+ if ti<0.0:
+ return 0.0
+ else:
+ curf=(maxf-minf)*ti/maxt+minf
+ if ti<(maxt*0.95):
+ return sin(2*pi*curf*ti)
+ else:
+ return 0.0
+
+
+def sweptsine(t,minf=0.0, maxf=10.0):
+ """Generate a sweptsine vector by calling Usweep for each ti in t."""
+ T=max(t)-min(t)
+ Us = [Usweep(ti,T,minf,maxf) for ti in t]
+ return array(Us)
+
+
+mytypes=['-','--',':','-.']
+colors=['b','y','r','g','c','k']#['y','b','r','g','c','k']
+
+def _getlinetype(ax=None):
+ if ax is None:
+ import pylab
+ ax = pylab.gca()
+ myind=ax._get_lines.count
+ return {'color':colors[myind % len(colors)],'linestyle':mytypes[myind % len(mytypes)]}
+
+
+def create_step_vector(t, step_time=0.0, amp=1.0):
+ u = zeros_like(t)
+ dt = t[1]-t[0]
+ indon = int(step_time/dt)
+ u[indon:] = amp
+ return u
+
+
+def rate_limiter(uin, du):
+ uout = zeros_like(uin)
+ N = len(uin)
+ for n in range(1,N):
+ curchange = uin[n]-uout[n-1]
+ if curchange > du:
+ uout[n] = uout[n-1]+du
+ elif curchange < -du:
+ uout[n] = uout[n-1]-du
+ else:
+ uout[n] = uin[n]
+ return uout
+
+
+
+
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: Ryan K. <rk...@si...> - 2010-06-17 21:19:05
|
That sounds great. It is good to know other people are working on this. I have had a few masters students ticker with it, but have been more or less been working on my own. I would be glad to be on both lists, but please use my list email: rya...@gm.... My TransferFunction class is basically a signal.ltisys with num and den as poly1d's (also from scipy). I mainly use the RootLocus and FreqResp methods along with the feedback function in my courses. > * You have implemented a couple of functions that I would love to pull into the python-control > library, including root locus plots and c2d support. Pulling those in would be a good way to see > how much synergy there is (based on how much rewrite is required). If this is OK by you, I > could give it a shot and let you know how it goes. The zoh c2d stuff took some work and I am happy with how it turned out (but I am not a digital control expert). So, this seems like a great first step. > * Going the other way, I've been trying to write the python-control code so that any object that > has the right attributes and/or member functions can be used. So, for example, the nyquist() >function should (eventually) work on any object for which there is a freqresp() member function > or a num and den attribute. I have tried to follow a similar approach where anything that has a num and den can be used in various parts of the code. For SISO systems, I think the code pretty much does what I would like it to. It could definitely use some well defined tests and better documentation and the scrutiny of more eyes on it, but I am more or less happy with it. The biggest need from my point of view is state-space stuff and especially MIMO. Actually, I don't even know if my code can elegantly handle SIMO. I need some SIMO tools for a project I am working on this summer at it is a little painful. I am certainly open to helping add new features. My rlocfind is slightly clumsy because I don't know who to click on a point on a plot and have IPython get the coordinates (without creating a wxPython gui just for the controls analysis or something). I will download your code in the next couple of days and give you my thoughts on how difficult it would be to combine the two code bases. -- Ryan Krauss, Ph.D. Assistant Professor Mechanical Engineering Southern Illinois University Edwardsville On Thu, Jun 17, 2010 at 3:43 PM, Richard Murray <mu...@cd...> wrote: > Hi Ryan. Good timing! I had just been trolling the web looking for others who are implementing control packages using python. Your name had already come up when I did a check a year or so back, but I was thinking that I would contact you to find out how things are going. And before I did, your message arrived! > > It would be great to find a way to combine efforts. I've got a couple of people who have expressed interest in adding to the functionality of the python-control package, so one possibility would be for us to take some of the functions that you have written and incorporate them into what we have so far. That would give us a sense of where the remaining gaps are and would allow us to look for some ways to continue to combine efforts. > > There are a couple of possible next steps: > > * You have implemented a couple of functions that I would love to pull into the python-control library, including root locus plots and c2d support. Pulling those in would be a good way to see how much synergy there is (based on how much rewrite is required). If this is OK by you, I could give it a shot and let you know how it goes. > > * Going the other way, I've been trying to write the python-control code so that any object that has the right attributes and/or member functions can be used. So, for example, the nyquist() function should (eventually) work on any object for which there is a freqresp() member function or a num and den attribute. It would be interesting to try this out and see whether we I can modify the python-control functions so that one of your TransferFunction objects could be passed to python-control functions. > > * Along these lines, if there are particular functions that you have been thinking about implementing next, please let us know so that we can either put them higher on the list or put them lower (if you will be working on them). We haven't yet done a prioritization, but my rough goal is to be able to run all of the examples in the controls book I wrote with Astrom (so that we can start using python in our classes as well), followed by some of the more advanced tools we use in a second quarter controls course (includes optimal control, trajectory generation, Kalman filter + extensions, etc). > > * Getting more active, if you would like to be signed up for either the python-control-announce list (low volume, summaries of when things are added) or the python-control-discuss (higher volume, used by people working on the project), just let me know. I'm cc'ing the latter list so that others know about the conversation. > > * Finally, if you actually have time to participate in the integration effort yourself, this would be wonderful. The "obvious" thing to do is to take the union of your TransferFunction class and the python-control TransferFunction class, then start working on generalizing various things to work for state space. For SISO this is trivial (signal.ltisys already keeps both num/den and ABCD representations for all objects), but gets tricker for things that might be MIMO. > > OK, that's probably enough for now -:). Let me know what you think about any and all of this. It would be great to collaborate! > > -richard > > On 17 Jun 2010, at 13:08 , Ryan Krauss wrote: > >> Richard, >> >> It looks like you are still actively developing a controls module for >> Python. I have used mine for teaching a first undergraduate controls >> course twice now and am relatively happy with it. It is poorly >> documented and doesn't really support state-space stuff yet, but it >> basically works for me. Let me know if you think it makes sense to >> combined our efforts somehow (I don't have a ton of time to spend on >> this, but would love to see other professors use Python for controls). >> >> Thanks, >> >> Ryan >> >> p.s. The version on my school webpage is slightly old. I need to fix >> that. Here are links to the temporary home of the latest versions >> from the course webpage: >> http://www.cs.siue.edu/~rkrauss/450/2010/controls_dist/controls-1.1.3.tar.gz >> http://www.cs.siue.edu/~rkrauss/450/2010/controls_dist/controls-1.1.3.win32.exe >> >> -- >> Ryan Krauss, Ph.D. >> Assistant Professor >> Mechanical Engineering >> Southern Illinois University Edwardsville > > |
|
From: Richard M. <mu...@cd...> - 2010-06-17 20:43:16
|
Hi Ryan. Good timing! I had just been trolling the web looking for others who are implementing control packages using python. Your name had already come up when I did a check a year or so back, but I was thinking that I would contact you to find out how things are going. And before I did, your message arrived! It would be great to find a way to combine efforts. I've got a couple of people who have expressed interest in adding to the functionality of the python-control package, so one possibility would be for us to take some of the functions that you have written and incorporate them into what we have so far. That would give us a sense of where the remaining gaps are and would allow us to look for some ways to continue to combine efforts. There are a couple of possible next steps: * You have implemented a couple of functions that I would love to pull into the python-control library, including root locus plots and c2d support. Pulling those in would be a good way to see how much synergy there is (based on how much rewrite is required). If this is OK by you, I could give it a shot and let you know how it goes. * Going the other way, I've been trying to write the python-control code so that any object that has the right attributes and/or member functions can be used. So, for example, the nyquist() function should (eventually) work on any object for which there is a freqresp() member function or a num and den attribute. It would be interesting to try this out and see whether we I can modify the python-control functions so that one of your TransferFunction objects could be passed to python-control functions. * Along these lines, if there are particular functions that you have been thinking about implementing next, please let us know so that we can either put them higher on the list or put them lower (if you will be working on them). We haven't yet done a prioritization, but my rough goal is to be able to run all of the examples in the controls book I wrote with Astrom (so that we can start using python in our classes as well), followed by some of the more advanced tools we use in a second quarter controls course (includes optimal control, trajectory generation, Kalman filter + extensions, etc). * Getting more active, if you would like to be signed up for either the python-control-announce list (low volume, summaries of when things are added) or the python-control-discuss (higher volume, used by people working on the project), just let me know. I'm cc'ing the latter list so that others know about the conversation. * Finally, if you actually have time to participate in the integration effort yourself, this would be wonderful. The "obvious" thing to do is to take the union of your TransferFunction class and the python-control TransferFunction class, then start working on generalizing various things to work for state space. For SISO this is trivial (signal.ltisys already keeps both num/den and ABCD representations for all objects), but gets tricker for things that might be MIMO. OK, that's probably enough for now -:). Let me know what you think about any and all of this. It would be great to collaborate! -richard On 17 Jun 2010, at 13:08 , Ryan Krauss wrote: > Richard, > > It looks like you are still actively developing a controls module for > Python. I have used mine for teaching a first undergraduate controls > course twice now and am relatively happy with it. It is poorly > documented and doesn't really support state-space stuff yet, but it > basically works for me. Let me know if you think it makes sense to > combined our efforts somehow (I don't have a ton of time to spend on > this, but would love to see other professors use Python for controls). > > Thanks, > > Ryan > > p.s. The version on my school webpage is slightly old. I need to fix > that. Here are links to the temporary home of the latest versions > from the course webpage: > http://www.cs.siue.edu/~rkrauss/450/2010/controls_dist/controls-1.1.3.tar.gz > http://www.cs.siue.edu/~rkrauss/450/2010/controls_dist/controls-1.1.3.win32.exe > > -- > Ryan Krauss, Ph.D. > Assistant Professor > Mechanical Engineering > Southern Illinois University Edwardsville |
|
From: <mur...@us...> - 2010-06-11 04:51:01
|
Revision: 24
http://python-control.svn.sourceforge.net/python-control/?rev=24&view=rev
Author: murrayrm
Date: 2010-06-11 04:50:54 +0000 (Fri, 11 Jun 2010)
Log Message:
-----------
minor code cleanup
Modified Paths:
--------------
trunk/examples/pvtol-lqr.py
trunk/examples/pvtol-nested.py
trunk/src/freqplot.py
Modified: trunk/examples/pvtol-lqr.py
===================================================================
--- trunk/examples/pvtol-lqr.py 2010-06-07 05:45:44 UTC (rev 23)
+++ trunk/examples/pvtol-lqr.py 2010-06-11 04:50:54 UTC (rev 24)
@@ -112,8 +112,6 @@
# Step response for the first input
H1ax = ss(Ax - Bx*K1a[0,lat], Bx*K1a[0,lat]*xd[lat,:], Cx, Dx);
(Tx, Yx) = step(H1ax, T=linspace(0,10,100));
-print Tx.shape
-print Yx.shape
# Step response for the second input
H1ay = ss(Ay - By*K1a[1,alt], By*K1a[1,alt]*yd[alt,:], Cy, Dy);
Modified: trunk/examples/pvtol-nested.py
===================================================================
--- trunk/examples/pvtol-nested.py 2010-06-07 05:45:44 UTC (rev 23)
+++ trunk/examples/pvtol-nested.py 2010-06-11 04:50:54 UTC (rev 24)
@@ -144,5 +144,5 @@
# Gang of Four
figure(11); clf();
-gangof4(Hi*Po, Co, linspace(-2, 3));
+gangof4(Hi*Po, Co);
Modified: trunk/src/freqplot.py
===================================================================
--- trunk/src/freqplot.py 2010-06-07 05:45:44 UTC (rev 23)
+++ trunk/src/freqplot.py 2010-06-11 04:50:54 UTC (rev 24)
@@ -52,7 +52,7 @@
Usage
=====
- (magh, phaseh) = bode(sys, omega=None, dB=False, Hz=False)
+ (magh, phaseh) = bode(syslist, omega=None, dB=False, Hz=False)
Plots a Bode plot for the system over a (optional) frequency range.
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mur...@us...> - 2010-06-07 05:45:51
|
Revision: 23
http://python-control.svn.sourceforge.net/python-control/?rev=23&view=rev
Author: murrayrm
Date: 2010-06-07 05:45:44 +0000 (Mon, 07 Jun 2010)
Log Message:
-----------
Working LQR example + lqr() fixes
* Fixed bugs in lqr related to matrix multiplication (* vs dot())
* Added PVTOL LQR example (based on AM08 example)
* Version 0.3c ready to be released...
Modified Paths:
--------------
trunk/ChangeLog
trunk/examples/pvtol-lqr.py
trunk/src/statefbk.py
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2010-06-06 21:24:07 UTC (rev 22)
+++ trunk/ChangeLog 2010-06-07 05:45:44 UTC (rev 23)
@@ -1,5 +1,7 @@
2010-06-06 Richard Murray <murray@sumatra.local>
+ * examples/pvtol-lqr.py: Added example to test out LQR routines
+
* src/matlab.py (bode): created a wrapper that allows MATLAB style
arguments for bode (eg, bode(sys1, sys2))
Modified: trunk/examples/pvtol-lqr.py
===================================================================
--- trunk/examples/pvtol-lqr.py 2010-06-06 21:24:07 UTC (rev 22)
+++ trunk/examples/pvtol-lqr.py 2010-06-07 05:45:44 UTC (rev 23)
@@ -1,192 +1,195 @@
-% pvtol_lqr.m - LQR design for vectored thrust aircraft
-% RMM, 14 Jan 03
+# pvtol_lqr.m - LQR design for vectored thrust aircraft
+# RMM, 14 Jan 03
+#
+# This file works through an LQR based design problem, using the
+# planar vertical takeoff and landing (PVTOL) aircraft example from
+# Astrom and Mruray, Chapter 5. It is intended to demonstrate the
+# basic functionality of the python-control package.
+#
-aminit;
+from numpy import * # Grab all of the NumPy functions
+from matplotlib.pyplot import * # Grab MATLAB plotting functions
+from control.matlab import * # MATLAB-like functions
-%%
-%% System dynamics
-%%
-%% These are the dynamics for the PVTOL system, written in state space
-%% form.
-%%
+#
+# System dynamics
+#
+# These are the dynamics for the PVTOL system, written in state space
+# form.
+#
-pvtol_params; % System parameters
+# System parameters
+m = 4; # mass of aircraft
+J = 0.0475; # inertia around pitch axis
+r = 0.25; # distance to center of force
+g = 9.8; # gravitational constant
+c = 0.05; # damping factor (estimated)
-% System matrices (entire plant: 2 input, 2 output)
-xe = [0 0 0 0 0 0]; ue = [0 m*g];
-[A, B, C, D] = pvtol_linearize(xe, ue);
+# State space dynamics
+xe = [0, 0, 0, 0, 0, 0]; # equilibrium point of interest
+ue = [0, m*g]; # (note these are lists, not matrices)
-%%
-%% Construct inputs and outputs corresponding to steps in xy position
-%%
-%% The vectors xd and yd correspond to the states that are the desired
-%% equilibrium states for the system. The matrices Cx and Cy are the
-%% corresponding outputs.
-%%
-%% The way these vectors are used is to compute the closed loop system
-%% dynamics as
-%%
-%% xdot = Ax + B u => xdot = (A-BK)x + K xd
-%% u = -K(x - xd) y = Cx
-%%
-%% The closed loop dynamics can be simulated using the "step" command,
-%% with K*xd as the input vector (assumes that the "input" is unit size,
-%% so that xd corresponds to the desired steady state.
-%%
+# Dynamics matrix (use matrix type so that * works for multiplication)
+A = matrix(
+ [[ 0, 0, 0, 1, 0, 0],
+ [ 0, 0, 0, 0, 1, 0],
+ [ 0, 0, 0, 0, 0, 1],
+ [ 0, 0, (-ue[0]*sin(xe[2]) - ue[1]*cos(xe[2]))/m, -c/m, 0, 0],
+ [ 0, 0, (ue[0]*cos(xe[2]) - ue[1]*sin(xe[2]))/m, 0, -c/m, 0],
+ [ 0, 0, 0, 0, 0, 0 ]])
-xd = [1; 0; 0; 0; 0; 0]; Cx = [1 0 0 0 0 0];
-yd = [0; 1; 0; 0; 0; 0]; Cy = [0 1 0 0 0 0];
+# Input matrix
+B = matrix(
+ [[0, 0], [0, 0], [0, 0],
+ [cos(xe[2])/m, -sin(xe[2])/m],
+ [sin(xe[2])/m, cos(xe[2])/m],
+ [r/J, 0]])
-%%
-%% LQR design
-%%
+# Output matrix
+C = matrix([[1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0]])
+D = matrix([[0, 0], [0, 0]])
-% Start with a diagonal weighting
-Qx1 = diag([1, 1, 1, 1, 1, 1]);
-Qu1a = diag([1, 1]);
-K1a = lqr(A, B, Qx1, Qu1a);
+#
+# Construct inputs and outputs corresponding to steps in xy position
+#
+# The vectors xd and yd correspond to the states that are the desired
+# equilibrium states for the system. The matrices Cx and Cy are the
+# corresponding outputs.
+#
+# The way these vectors are used is to compute the closed loop system
+# dynamics as
+#
+# xdot = Ax + B u => xdot = (A-BK)x + K xd
+# u = -K(x - xd) y = Cx
+#
+# The closed loop dynamics can be simulated using the "step" command,
+# with K*xd as the input vector (assumes that the "input" is unit size,
+# so that xd corresponds to the desired steady state.
+#
-% Close the loop: xdot = Ax + B K (x-xd)
-H1a = ss(A-B*K1a, B*K1a*[xd, yd], [Cx; Cy], 0);
-[Y, T] = step(H1a, 10);
+xd = matrix([[1], [0], [0], [0], [0], [0]]);
+yd = matrix([[0], [1], [0], [0], [0], [0]]);
-figure(1); clf; subplot(321);
- plot(T, Y(:,1, 1), '-', T, Y(:,2, 2), '--', ...
- 'Linewidth', AM_data_linewidth); hold on;
- plot([0 10], [1 1], 'k-', 'Linewidth', AM_ref_linewidth); hold on;
+#
+# Extract the relevant dynamics for use with SISO library
+#
+# The current python-control library only supports SISO transfer
+# functions, so we have to modify some parts of the original MATLAB
+# code to extract out SISO systems. To do this, we define the 'lat' and
+# 'alt' index vectors to consist of the states that are are relevant
+# to the lateral (x) and vertical (y) dynamics.
+#
- amaxis([0, 10, -0.1, 1.4]);
- xlabel('time'); ylabel('position');
+# Indices for the parts of the state that we want
+lat = (0,2,3,5);
+alt = (1,4);
- lgh = legend('x', 'y', 'Location', 'southeast');
- legend(lgh, 'boxoff');
- amprint('pvtol-lqrstep1.eps');
+# Decoupled dynamics
+Ax = (A[lat, :])[:, lat]; #! not sure why I have to do it this way
+Bx = B[lat, 0]; Cx = C[0, lat]; Dx = D[0, 0];
-% Look at different input weightings
-Qu1a = diag([1, 1]); K1a = lqr(A, B, Qx1, Qu1a);
-H1ax = ss(A-B*K1a,B(:,1)*K1a(1,:)*xd,Cx,0);
+Ay = (A[alt, :])[:, alt]; #! not sure why I have to do it this way
+By = B[alt, 1]; Cy = C[1, alt]; Dy = D[1, 1];
-Qu1b = 40^2*diag([1, 1]); K1b = lqr(A, B, Qx1, Qu1b);
-H1bx = ss(A-B*K1b,B(:,1)*K1b(1,:)*xd,Cx,0);
+# Label the plot
+clf();
+suptitle("LQR controllers for vectored thrust aircraft (pvtol-lqr)")
-Qu1c = 200^2*diag([1, 1]); K1c = lqr(A, B, Qx1, Qu1c);
-H1cx = ss(A-B*K1c,B(:,1)*K1c(1,:)*xd,Cx,0);
+#
+# LQR design
+#
-[Y1, T1] = step(H1ax, 10);
-[Y2, T2] = step(H1bx, 10);
-[Y3, T3] = step(H1cx, 10);
+# Start with a diagonal weighting
+Qx1 = diag([1, 1, 1, 1, 1, 1]);
+Qu1a = diag([1, 1]);
+(K, X, E) = lqr(A, B, Qx1, Qu1a); K1a = matrix(K);
-figure(2); clf; subplot(321);
- plot(T1, Y1, 'b-', 'Linewidth', AM_data_linewidth); hold on;
- plot(T2, Y2, 'b-', 'Linewidth', AM_data_linewidth); hold on;
- plot(T3, Y3, 'b-', 'Linewidth', AM_data_linewidth); hold on;
- plot([0 10], [1 1], 'k-', 'Linewidth', AM_ref_linewidth); hold on;
+# Close the loop: xdot = Ax - B K (x-xd)
+# Note: python-control requires we do this 1 input at a time
+# H1a = ss(A-B*K1a, B*K1a*concatenate((xd, yd), axis=1), C, D);
+# (T, Y) = step(H1a, T=linspace(0,10,100));
- amaxis([0, 10, -0.1, 1.4]);
- xlabel('time'); ylabel('position');
+# Step response for the first input
+H1ax = ss(Ax - Bx*K1a[0,lat], Bx*K1a[0,lat]*xd[lat,:], Cx, Dx);
+(Tx, Yx) = step(H1ax, T=linspace(0,10,100));
+print Tx.shape
+print Yx.shape
- arcarrow([1.3 0.8], [5 0.45], -6);
- text(5.3, 0.4, 'rho');
+# Step response for the second input
+H1ay = ss(Ay - By*K1a[1,alt], By*K1a[1,alt]*yd[alt,:], Cy, Dy);
+(Ty, Yy) = step(H1ay, T=linspace(0,10,100));
- amprint('pvtol-lqrstep2.eps');
+subplot(221); title("Identity weights")
+# plot(T, Y[:,1, 1], '-', T, Y[:,2, 2], '--'); hold(True);
+plot(Tx, Yx[0,:].T, '-', Ty, Yy[0,:].T, '--'); hold(True);
+plot([0, 10], [1, 1], 'k-'); hold(True);
-% Output weighting - change Qx to use outputs
-Qx2 = [Cx; Cy]' * [Cx; Cy];
-Qu2 = 0.1 * diag([1, 1]);
-K2 = lqr(A, B, Qx2, Qu2);
+axis([0, 10, -0.1, 1.4]);
+ylabel('position');
+legend(('x', 'y'), loc='lower right');
-H2x = ss(A-B*K2,B(:,1)*K2(1,:)*xd,Cx,0);
-H2y = ss(A-B*K2,B(:,2)*K2(2,:)*yd,Cy,0);
+# Look at different input weightings
+Qu1a = diag([1, 1]); (K1a, X, E) = lqr(A, B, Qx1, Qu1a);
+H1ax = ss(Ax - Bx*K1a[0,lat], Bx*K1a[0,lat]*xd[lat,:], Cx, Dx);
-figure(3); step(H2x, H2y, 10);
-legend('x', 'y');
+Qu1b = (40**2)*diag([1, 1]); (K1b, X, E) = lqr(A, B, Qx1, Qu1b);
+H1bx = ss(Ax - Bx*K1b[0,lat], Bx*K1b[0,lat]*xd[lat,:],Cx, Dx);
-%%
-%% Physically motivated weighting
-%%
-%% Shoot for 1 cm error in x, 10 cm error in y. Try to keep the angle
-%% less than 5 degrees in making the adjustments. Penalize side forces
-%% due to loss in efficiency.
-%%
+Qu1c = (200**2)*diag([1, 1]); (K1c, X, E) = lqr(A, B, Qx1, Qu1c);
+H1cx = ss(Ax - Bx*K1c[0,lat], Bx*K1c[0,lat]*xd[lat,:],Cx, Dx);
-Qx3 = diag([100, 10, 2*pi/5, 0, 0, 0]);
-Qu3 = 0.1 * diag([1, 10]);
-K3 = lqr(A, B, Qx3, Qu3);
+[T1, Y1] = step(H1ax, T=linspace(0,10,100));
+[T2, Y2] = step(H1bx, T=linspace(0,10,100));
+[T3, Y3] = step(H1cx, T=linspace(0,10,100));
-H3x = ss(A-B*K3,B(:,1)*K3(1,:)*xd,Cx,0);
-H3y = ss(A-B*K3,B(:,2)*K3(2,:)*yd,Cy,0);
-figure(4); clf; subplot(221);
-step(H3x, H3y, 10);
-legend('x', 'y');
+subplot(222); title("Effect of input weights")
+plot(T1, Y1[0,:].T, 'b-'); hold(True);
+plot(T2, Y2[0,:].T, 'b-'); hold(True);
+plot(T3, Y3[0,:].T, 'b-'); hold(True);
+plot([0 ,10], [1, 1], 'k-'); hold(True);
-%%
-%% Velocity control
-%%
-%% In this example, we modify the system so that we control the
-%% velocity of the system in the x direction. We ignore the
-%% dynamics in the vertical (y) direction. These dynamics demonstrate
-%% the role of the feedforward system since the equilibrium point
-%% corresponding to vd neq 0 requires a nonzero input.
-%%
-%% For this example, we use a control law u = -K(x-xd) + ud and convert
-%% this to the form u = -K x + N r, where r is the reference input and
-%% N is computed as described in class.
-%%
+axis([0, 10, -0.1, 1.4]);
-% Extract system dynamics: theta, xdot, thdot
-Av = A([3 4 6], [3 4 6]);
-Bv = B([3 4 6], 1);
-Cv = [0 1 0]; % choose vx as output
-Dv = 0;
+# arcarrow([1.3, 0.8], [5, 0.45], -6);
+text(5.3, 0.4, 'rho');
-% Design the feedback term using LQR
-Qxv = diag([2*pi/5, 10, 0]);
-Quv = 0.1;
-Kv = lqr(Av, Bv, Qxv, Quv);
+# Output weighting - change Qx to use outputs
+Qx2 = C.T * C;
+Qu2 = 0.1 * diag([1, 1]);
+(K, X, E) = lqr(A, B, Qx2, Qu2); K2 = matrix(K)
-% Design the feedforward term by solve for eq pt in terms of reference r
-T = [Av Bv; Cv Dv]; % system matrix
-Nxu = T \ [0; 0; 0; 1]; % compute [Nx; Nu]
-Nx = Nxu(1:3); Nu = Nxu(4); % extract Nx and Nu
-N = Nu + Kv*Nx; % compute feedforward term
+H2x = ss(Ax - Bx*K2[0,lat], Bx*K2[0,lat]*xd[lat,:], Cx, Dx);
+H2y = ss(Ay - By*K2[1,alt], By*K2[1,alt]*yd[alt,:], Cy, Dy);
-%%
-%% Design #1: no feedforward input, ud
-%%
+subplot(223); title("Output weighting")
+[T2x, Y2x] = step(H2x, T=linspace(0,10,100));
+[T2y, Y2y] = step(H2y, T=linspace(0,10,100));
+plot(T2x, Y2x[0,:].T, T2y, Y2y[0,:].T)
+ylabel('position');
+xlabel('time'); ylabel('position');
+legend(('x', 'y'), loc='lower right');
-Nv1 = [0; 1; 0];
-Hv1 = ss(Av-Bv*Kv, Bv*Kv*Nx, Cv, 0);
-step(Hv1, 10);
+#
+# Physically motivated weighting
+#
+# Shoot for 1 cm error in x, 10 cm error in y. Try to keep the angle
+# less than 5 degrees in making the adjustments. Penalize side forces
+# due to loss in efficiency.
+#
-%%
-%% Design #2: compute feedforward gain corresponding to equilibrium point
-%%
+Qx3 = diag([100, 10, 2*pi/5, 0, 0, 0]);
+Qu3 = 0.1 * diag([1, 10]);
+(K, X, E) = lqr(A, B, Qx3, Qu3); K3 = matrix(K);
-Hv2 = ss(Av-Bv*Kv, Bv*N, Cv, 0);
-step(Hv2, 10);
+H3x = ss(Ax - Bx*K3[0,lat], Bx*K3[0,lat]*xd[lat,:], Cx, Dx);
+H3y = ss(Ay - By*K3[1,alt], By*K3[1,alt]*yd[alt,:], Cy, Dy);
+subplot(224)
+# step(H3x, H3y, 10);
+[T3x, Y3x] = step(H3x, T=linspace(0,10,100));
+[T3y, Y3y] = step(H3y, T=linspace(0,10,100));
+plot(T3x, Y3x[0,:].T, T3y, Y3y[0,:].T)
+title("Physically motivated weights")
+xlabel('time');
+legend(('x', 'y'), loc='lower right');
-%%
-%% Design #3: integral action
-%%
-%% Add a new state to the system that is given by xidot = v - vd. We
-%% construct the control law by computing an LQR gain for the augmented
-%% system.
-%%
-
-Ai = [Av, [0; 0; 0]; [Cv, 0]];
-Bi = [Bv; 0];
-Ci = [Cv, 0];
-Di = Dv;
-
-% Design the feedback term, including weight on integrator error
-Qxi = diag([2*pi/5, 10, 0, 10]);
-Qui = 0.1;
-Ki = lqr(Ai, Bi, Qxi, Qui);
-
-% Desired state (augmented)
-xid = [0; 1; 0; 0];
-
-% Construct the closed loop system (including integrator)
-Hi = ss(Ai-Bi*Ki,Bi*Ki*xid - [0; 0; 0; Ci*xid],Ci,0);
-step(Hi, 10);
-
+show()
Modified: trunk/src/statefbk.py
===================================================================
--- trunk/src/statefbk.py 2010-06-06 21:24:07 UTC (rev 22)
+++ trunk/src/statefbk.py 2010-06-07 05:45:44 UTC (rev 23)
@@ -41,6 +41,7 @@
# External packages and modules
import numpy as np
+import ctrlutil
from control.exception import *
# Pole placement
@@ -139,8 +140,7 @@
if (len(args) < 4):
raise ControlArgument("not enough input arguments")
- elif (getattr(args[0], 'A', None) and
- getattr(args[0], 'B', None)):
+ elif (ctrlutil.issys(args[0])):
# We were passed a system as the first argument; extract A and B
#! TODO: really just need to check for A and B attributes
A = np.array(args[0].A, ndmin=2, dtype=float);
@@ -179,7 +179,7 @@
X,rcond,w,S,U = sb02md(nstates, A_b, G, Q_b, 'C')
# Now compute the return value
- K = np.linalg.inv(R) * (np.transpose(B) * X + np.transpose(N));
+ K = np.dot(np.linalg.inv(R), (np.dot(B.T, X) + N.T));
S = X;
E = w[0:nstates];
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: Tony V. <ton...@al...> - 2010-06-06 21:45:57
|
Got it this time. -tv On Jun 6, 2010, at 14:25, Richard Murray <mu...@cd...> wrote: > OK, I've now set it up so that SVN commit messages go to the discuss > list, so we should all see them. > > -richard > > On 6 Jun 2010, at 14:11 , Tony Vanelli wrote: > >> Doesn't look like I received it. I also checked my spam folder just >> in >> case. >> >> -tv >> >> >> On Jun 6, 2010, at 10:31, Richard Murray <mu...@cd...> >> wrote: >> >>> Tony and Erik: >>> >>> Can you let me know if you got the note below from SourceForge. I >>> just posted some changes to the package (needed them for an example >>> I'm putting together) and want to make sure that everyone is getting >>> copies of the changes that got posted. If not, I'll sort out what >>> we need to change to make sure you get these sorts of e-mails. >>> >>> Thanks. >>> >>> -richard >>> >>> Begin forwarded message: >>> >>>> From: mur...@us... >>>> Date: 6 June 2010 10:29:43 PDT >>>> To: mu...@cd... >>>> Subject: SF.net SVN: python-control:[21] trunk >>>> >>>> Revision: 21 >>>> http://python-control.svn.sourceforge.net/python-control/?rev=21&view=rev >>>> Author: murrayrm >>>> Date: 2010-06-06 17:29:42 +0000 (Sun, 06 Jun 2010) >>>> >>>> Log Message: >>>> ----------- >>>> Added changes to allow bode() to plot frequency resp for multiple >>>> systems >>>> * freqplot.bode now allows a list as the first argument >>>> * Added matlab.bode wrapper for MATLAB syntax: bode(sys1, ..., >>>> sysN) >>>> >>>> Added automatic frequency range determination in bode >>>> * Computes range of frequencies based on pole/zero locations >>>> * Eventually need to split this out as a ctrlutil function for use >>>> elsewhere >>>> >>>> Modified Paths: >>>> -------------- >>>> trunk/ChangeLog >>>> trunk/src/ctrlutil.py >>>> trunk/src/freqplot.py >>>> trunk/src/matlab.py >>>> trunk/src/pzmap.py >>>> trunk/src/statesp.py >>>> trunk/src/xferfcn.py >>>> >>>> Modified: trunk/ChangeLog >>>> =================================================================== >>>> --- trunk/ChangeLog 2010-06-01 02:08:20 UTC (rev 20) >>>> +++ trunk/ChangeLog 2010-06-06 17:29:42 UTC (rev 21) >>>> @@ -1,3 +1,25 @@ >>>> +2010-06-06 Richard Murray <murray@sumatra.local> >>>> + >>>> + * src/matlab.py (bode): created a wrapper that allows MATLAB >>>> style >>>> + arguments for bode (eg, bode(sys1, sys2)) >>>> + >>>> + * src/ctrlutil.py (issys): added new function to check if an >>>> object >>>> + is a system (state space or transfer function). Will >>>> generalize >>>> + this latter to look for other compatible classes. >>>> + >>>> + * src/freqplot.py (bode): Compute frequency range of bode plot >>>> based >>>> + on poles and zeros >>>> + (bode): Allow bode plot to be passed a list (or tuple) as the >>>> first >>>> + argument, in which case multiple bode plots are generated >>>> + >>>> + * src/statesp.py (StateSpace.zeros): new function to compute >>>> zeros >>>> + for a state space system >>>> + (StateSpace): defined new functions to compute poles of a >>>> state >>>> + space system >>>> + >>>> + * src/xferfcn.py (TransferFunction): defined new functions to >>>> + compute poles and zeros of a transfer function. >>>> + >>>> 2010-05-31 Richard Murray <murray@sumatra.local> >>>> >>>> * src/exception.py (ControlNotImplemented): added new exception, >>>> to >>>> >>>> Modified: trunk/src/ctrlutil.py >>>> =================================================================== >>>> --- trunk/src/ctrlutil.py 2010-06-01 02:08:20 UTC (rev 20) >>>> +++ trunk/src/ctrlutil.py 2010-06-06 17:29:42 UTC (rev 21) >>>> @@ -40,6 +40,12 @@ >>>> # >>>> # $Id$ >>>> >>>> +# Packages that we need access to >>>> +import scipy as sp >>>> +import statesp >>>> +import xferfcn >>>> + >>>> +# Specific functions that we use >>>> from scipy import pi >>>> >>>> # Utility function to unwrap an angle measurement >>>> @@ -78,3 +84,12 @@ >>>> >>>> # return the updated list >>>> return angle >>>> + >>>> +# Determine if an object is a system >>>> +def issys(object): >>>> + # Check for a member of one of the classes that we define here >>>> + if (isinstance(object, (statesp.StateSpace, >>>> xferfcn.TransferFunction))): >>>> + return True >>>> + >>>> + # Didn't find anything that matched >>>> + return False >>>> >>>> Modified: trunk/src/freqplot.py >>>> =================================================================== >>>> --- trunk/src/freqplot.py 2010-06-01 02:08:20 UTC (rev 20) >>>> +++ trunk/src/freqplot.py 2010-06-06 17:29:42 UTC (rev 21) >>>> @@ -42,11 +42,12 @@ >>>> >>>> import matplotlib.pyplot as plt >>>> import scipy as sp >>>> +import numpy as np >>>> from ctrlutil import unwrap >>>> from bdalg import feedback >>>> >>>> # Bode plot >>>> -def bode(sys, omega=None, dB=False, Hz=False): >>>> +def bode(syslist, omega=None, dB=False, Hz=False): >>>> """Bode plot for a system >>>> >>>> Usage >>>> @@ -57,8 +58,8 @@ >>>> >>>> Parameters >>>> ---------- >>>> - sys : linsys >>>> - Linear input/output system >>>> + syslist : linsys >>>> + List of linear input/output systems (single system is OK) >>>> omega : freq_range >>>> Range of frequencies (list or bounds) in rad/sec >>>> dB : boolean >>>> @@ -76,45 +77,80 @@ >>>> 1. Use (mag, phase, freq) = sys.freqresp(freq) to generate the >>>> frequency response for a system. >>>> """ >>>> + # If argument was a singleton, turn it into a list >>>> + if (not getattr(syslist, '__iter__', False)): >>>> + syslist = (syslist,) >>>> + >>>> + # >>>> # Select a default range if none is provided >>>> + # >>>> + # This code looks at the poles and zeros of all of the systems >>>> that >>>> + # we are plotting and sets the frequency range to be one >>>> decade above >>>> + # and below the min and max feature frequencies, rounded to >>>> the nearest >>>> + # integer. It excludes poles and zeros at the origin. If no >>>> features >>>> + # are found, it turns logspace(-1, 1) >>>> + # >>>> if (omega == None): >>>> - omega = sp.logspace(-2, 2); >>>> + # Find the list of all poles and zeros in the systems >>>> + features = np.array(()) >>>> + for sys in syslist: >>>> + # Add new features to the list >>>> + features = np.concatenate((features, np.abs >>>> (sys.poles))) >>>> + features = np.concatenate((features, np.abs >>>> (sys.zeros))) >>>> >>>> - # Get the magnitude and phase of the system >>>> - mag, phase, omega = sys.freqresp(omega) >>>> - if Hz: >>>> - omega = omega/(2*sp.pi) >>>> - if dB: >>>> - mag = 20*sp.log10(mag) >>>> - phase = unwrap(phase*180/sp.pi, 360) >>>> + # Get rid of poles and zeros at the origin >>>> + features = features[features != 0]; >>>> >>>> - # Get the dimensions of the current axis, which we will divide >>>> up >>>> - #! TODO: Not current implemented; just use subplot for now >>>> + # Make sure there is at least one point in the range >>>> + if (features.shape[0] == 0): features = [1]; >>>> >>>> - # Magnitude plot >>>> - plt.subplot(211); >>>> - if dB: >>>> - plt.semilogx(omega, mag) >>>> - else: >>>> - plt.loglog(omega, mag) >>>> - plt.grid(True) >>>> - plt.grid(True, which='minor') >>>> - if dB: >>>> - plt.ylabel("Magnitude (dB)") >>>> - else: >>>> - plt.ylabel("Magnitude") >>>> + # Take the log of the features >>>> + features = np.log10(features) >>>> >>>> - # Phase plot >>>> - plt.subplot(212); >>>> - plt.semilogx(omega, phase) >>>> - plt.grid(True) >>>> - plt.grid(True, which='minor') >>>> - plt.ylabel("Phase (deg)") >>>> - if Hz: >>>> - plt.xlabel("Frequency (Hz)") >>>> - else: >>>> - plt.xlabel("Frequency (rad/sec)") >>>> + # Set the to be an order of magnitude beyond any features >>>> + omega = sp.logspace(np.floor(np.min(features))-1, >>>> + np.ceil(np.max(features))+1) >>>> >>>> + for sys in syslist: >>>> + # Get the magnitude and phase of the system >>>> + mag, phase, omega = sys.freqresp(omega) >>>> + if Hz: omega = omega/(2*sp.pi) >>>> + if dB: mag = 20*sp.log10(mag) >>>> + phase = unwrap(phase*180/sp.pi, 360) >>>> + >>>> + # Get the dimensions of the current axis, which we will >>>> divide up >>>> + #! TODO: Not current implemented; just use subplot for now >>>> + >>>> + # Magnitude plot >>>> + plt.subplot(211); >>>> + if dB: >>>> + plt.semilogx(omega, mag) >>>> + plt.ylabel("Magnitude (dB)") >>>> + else: >>>> + plt.loglog(omega, mag) >>>> + plt.ylabel("Magnitude") >>>> + >>>> + # Add a grid to the plot >>>> + plt.grid(True) >>>> + plt.grid(True, which='minor') >>>> + plt.hold(True); >>>> + >>>> + # Phase plot >>>> + plt.subplot(212); >>>> + plt.semilogx(omega, phase) >>>> + plt.hold(True) >>>> + >>>> + # Add a grid to the plot >>>> + plt.grid(True) >>>> + plt.grid(True, which='minor') >>>> + plt.ylabel("Phase (deg)") >>>> + >>>> + # Label the frequency axis >>>> + if Hz: >>>> + plt.xlabel("Frequency (Hz)") >>>> + else: >>>> + plt.xlabel("Frequency (rad/sec)") >>>> + >>>> return (211, 212) >>>> >>>> # Nyquist plot >>>> >>>> Modified: trunk/src/matlab.py >>>> =================================================================== >>>> --- trunk/src/matlab.py 2010-06-01 02:08:20 UTC (rev 20) >>>> +++ trunk/src/matlab.py 2010-06-06 17:29:42 UTC (rev 21) >>>> @@ -50,6 +50,7 @@ >>>> >>>> # Libraries that we make use of >>>> import scipy as sp # SciPy library (used all over) >>>> +import numpy as np # NumPy library >>>> import scipy.signal as signal # Signal processing library >>>> >>>> # Import MATLAB-like functions that are defined in other packages >>>> @@ -57,11 +58,16 @@ >>>> from scipy.signal import lsim, impulse, step >>>> from scipy import linspace, logspace >>>> >>>> -# Import MATLAB-like functions that belong elsewhere in python- >>>> control >>>> -from ctrlutil import unwrap >>>> -from freqplot import bode, nyquist, gangof4 >>>> +# Control system library >>>> +import ctrlutil >>>> +import freqplot >>>> from statesp import StateSpace >>>> from xferfcn import TransferFunction >>>> +from exception import * >>>> + >>>> +# Import MATLAB-like functions that can be used as-is >>>> +from ctrlutil import unwrap >>>> +from freqplot import nyquist, gangof4 >>>> from bdalg import series, parallel, negate, feedback >>>> from pzmap import pzmap >>>> from statefbk import place, lqr >>>> @@ -304,3 +310,63 @@ >>>> def freqresp(H, omega): >>>> """Return the frequency response for an object H at frequency >>>> omega""" >>>> return H.freqresp(omega) >>>> + >>>> +# Bode plots >>>> +def bode(*args, **keywords): >>>> + """Bode plot of the frequency response >>>> + >>>> + Usage >>>> + ===== >>>> + bode(sys) >>>> + bode(sys, w) >>>> + bode(sys1, sys2, ..., sysN) >>>> + bode(sys1, sys2, ..., sysN, w) >>>> + bode(sys1, 'plotstyle1', ..., sysN, 'plotstyleN') >>>> + """ >>>> + >>>> + # If the first argument is a list, then assume python-control >>>> calling format >>>> + if (getattr(args[0], '__iter__', False)): >>>> + return freqplot.bode(*args, **keywords) >>>> + >>>> + # Otherwise, run through the arguments and collect up >>>> arguments >>>> + syslist = []; plotstyle=[]; omega=None; >>>> + i = 0; >>>> + while i < len(args): >>>> + # Check to see if this is a system of some sort >>>> + if (ctrlutil.issys(args[i])): >>>> + # Append the system to our list of systems >>>> + syslist.append(args[i]) >>>> + i += 1 >>>> + >>>> + # See if the next object is a plotsytle (string) >>>> + if (i < len(args) and isinstance(args[i], str)): >>>> + plotstyle.append(args[i]) >>>> + i += 1 >>>> + >>>> + # Go on to the next argument >>>> + continue >>>> + >>>> + # See if this is a frequency list >>>> + elif (isinstance(args[i], (list, np.ndarray))): >>>> + omega = args[i] >>>> + i += 1 >>>> + break >>>> + >>>> + else: >>>> + raise ControlArgument("unrecognized argument type") >>>> + >>>> + # Check to make sure that we processed all arguments >>>> + if (i < len(args)): >>>> + raise ControlArgument("not all arguments processed") >>>> + >>>> + # Check to make sure we got the same number of plotstyles as >>>> systems >>>> + if (len(plotstyle) != 0 and len(syslist) != len(plotstyle)): >>>> + raise ControlArgument("number of systems and plotstyles >>>> should be equal") >>>> + >>>> + # Warn about unimplemented plotstyles >>>> + #! TODO: remove this when plot styles are implemented in bode >>>> () >>>> + if (len(plotstyle) != 0): >>>> + print("Warning (matabl.bode): plot styles not >>>> implemented"); >>>> + >>>> + # Call the bode command >>>> + return freqplot.bode(syslist, omega, **keywords) >>>> >>>> Modified: trunk/src/pzmap.py >>>> =================================================================== >>>> --- trunk/src/pzmap.py 2010-06-01 02:08:20 UTC (rev 20) >>>> +++ trunk/src/pzmap.py 2010-06-06 17:29:42 UTC (rev 21) >>>> @@ -45,6 +45,7 @@ >>>> import xferfcn >>>> >>>> # Compute poles and zeros for a system >>>> +#! TODO: extend this to handle state space systems >>>> def pzmap(sys, Plot=True): >>>> """Plot a pole/zero map for a transfer function""" >>>> if (isinstance(sys, xferfcn.TransferFunction)): >>>> >>>> Modified: trunk/src/statesp.py >>>> =================================================================== >>>> --- trunk/src/statesp.py 2010-06-01 02:08:20 UTC (rev 20) >>>> +++ trunk/src/statesp.py 2010-06-06 17:29:42 UTC (rev 21) >>>> @@ -98,6 +98,18 @@ >>>> #! TODO: Not implemented >>>> return None >>>> >>>> + # Compute poles and zeros >>>> + def poles(self): return sp.roots(sp.poly(self.A)) >>>> + def zeros(self): >>>> + den = sp.poly1d(sp.poly(self.A)) >>>> + >>>> + # Compute the numerator based on zeros >>>> + #! TODO: This is currently limited to SISO systems >>>> + num = sp.poly1d(\ >>>> + sp.poly(self.A - sp.dot(self.B, self.C)) + (self.D[0] >>>> - 1) * den) >>>> + >>>> + return (sp.roots(num)) >>>> + >>>> # Negation of a system >>>> def __neg__(self): >>>> """Negate a state space system""" >>>> >>>> Modified: trunk/src/xferfcn.py >>>> =================================================================== >>>> --- trunk/src/xferfcn.py 2010-06-01 02:08:20 UTC (rev 20) >>>> +++ trunk/src/xferfcn.py 2010-06-06 17:29:42 UTC (rev 21) >>>> @@ -190,6 +190,10 @@ >>>> >>>> return mag, phase, omega >>>> >>>> + # Compute poles and zeros >>>> + def poles(self): return sp.roots(self.den) >>>> + def zeros(self): return sp.roots(self.num) >>>> + >>>> # Feedback around a transfer function >>>> def feedback(sys1, sys2, sign=-1): >>>> """Feedback interconnection between two transfer functions""" >>>> >>>> >>>> This was sent by the SourceForge.net collaborative development >>>> platform, the world's largest Open Source development site. >>> >>> >>> --- >>> --- >>> --- >>> --- >>> ------------------------------------------------------------------ >>> ThinkGeek and WIRED's GeekDad team up for the Ultimate >>> GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the >>> lucky parental unit. See the prize list and enter to win: >>> http://p.sf.net/sfu/thinkgeek-promo >>> _______________________________________________ >>> python-control-discuss mailing list >>> pyt...@li... >>> https://lists.sourceforge.net/lists/listinfo/python-control-discuss >> >> --- >> --- >> --- >> --------------------------------------------------------------------- >> ThinkGeek and WIRED's GeekDad team up for the Ultimate >> GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the >> lucky parental unit. See the prize list and enter to win: >> http://p.sf.net/sfu/thinkgeek-promo >> _______________________________________________ >> python-control-discuss mailing list >> pyt...@li... >> https://lists.sourceforge.net/lists/listinfo/python-control-discuss > > > --- > --- > --- > --------------------------------------------------------------------- > ThinkGeek and WIRED's GeekDad team up for the Ultimate > GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the > lucky parental unit. See the prize list and enter to win: > http://p.sf.net/sfu/thinkgeek-promo > _______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss |
|
From: Richard M. <mu...@cd...> - 2010-06-06 21:25:51
|
OK, I've now set it up so that SVN commit messages go to the discuss list, so we should all see them. -richard On 6 Jun 2010, at 14:11 , Tony Vanelli wrote: > Doesn't look like I received it. I also checked my spam folder just in > case. > > -tv > > > On Jun 6, 2010, at 10:31, Richard Murray <mu...@cd...> wrote: > >> Tony and Erik: >> >> Can you let me know if you got the note below from SourceForge. I >> just posted some changes to the package (needed them for an example >> I'm putting together) and want to make sure that everyone is getting >> copies of the changes that got posted. If not, I'll sort out what >> we need to change to make sure you get these sorts of e-mails. >> >> Thanks. >> >> -richard >> >> Begin forwarded message: >> >>> From: mur...@us... >>> Date: 6 June 2010 10:29:43 PDT >>> To: mu...@cd... >>> Subject: SF.net SVN: python-control:[21] trunk >>> >>> Revision: 21 >>> http://python-control.svn.sourceforge.net/python-control/?rev=21&view=rev >>> Author: murrayrm >>> Date: 2010-06-06 17:29:42 +0000 (Sun, 06 Jun 2010) >>> >>> Log Message: >>> ----------- >>> Added changes to allow bode() to plot frequency resp for multiple >>> systems >>> * freqplot.bode now allows a list as the first argument >>> * Added matlab.bode wrapper for MATLAB syntax: bode(sys1, ..., sysN) >>> >>> Added automatic frequency range determination in bode >>> * Computes range of frequencies based on pole/zero locations >>> * Eventually need to split this out as a ctrlutil function for use >>> elsewhere >>> >>> Modified Paths: >>> -------------- >>> trunk/ChangeLog >>> trunk/src/ctrlutil.py >>> trunk/src/freqplot.py >>> trunk/src/matlab.py >>> trunk/src/pzmap.py >>> trunk/src/statesp.py >>> trunk/src/xferfcn.py >>> >>> Modified: trunk/ChangeLog >>> =================================================================== >>> --- trunk/ChangeLog 2010-06-01 02:08:20 UTC (rev 20) >>> +++ trunk/ChangeLog 2010-06-06 17:29:42 UTC (rev 21) >>> @@ -1,3 +1,25 @@ >>> +2010-06-06 Richard Murray <murray@sumatra.local> >>> + >>> + * src/matlab.py (bode): created a wrapper that allows MATLAB >>> style >>> + arguments for bode (eg, bode(sys1, sys2)) >>> + >>> + * src/ctrlutil.py (issys): added new function to check if an >>> object >>> + is a system (state space or transfer function). Will generalize >>> + this latter to look for other compatible classes. >>> + >>> + * src/freqplot.py (bode): Compute frequency range of bode plot >>> based >>> + on poles and zeros >>> + (bode): Allow bode plot to be passed a list (or tuple) as the >>> first >>> + argument, in which case multiple bode plots are generated >>> + >>> + * src/statesp.py (StateSpace.zeros): new function to compute >>> zeros >>> + for a state space system >>> + (StateSpace): defined new functions to compute poles of a state >>> + space system >>> + >>> + * src/xferfcn.py (TransferFunction): defined new functions to >>> + compute poles and zeros of a transfer function. >>> + >>> 2010-05-31 Richard Murray <murray@sumatra.local> >>> >>> * src/exception.py (ControlNotImplemented): added new exception, >>> to >>> >>> Modified: trunk/src/ctrlutil.py >>> =================================================================== >>> --- trunk/src/ctrlutil.py 2010-06-01 02:08:20 UTC (rev 20) >>> +++ trunk/src/ctrlutil.py 2010-06-06 17:29:42 UTC (rev 21) >>> @@ -40,6 +40,12 @@ >>> # >>> # $Id$ >>> >>> +# Packages that we need access to >>> +import scipy as sp >>> +import statesp >>> +import xferfcn >>> + >>> +# Specific functions that we use >>> from scipy import pi >>> >>> # Utility function to unwrap an angle measurement >>> @@ -78,3 +84,12 @@ >>> >>> # return the updated list >>> return angle >>> + >>> +# Determine if an object is a system >>> +def issys(object): >>> + # Check for a member of one of the classes that we define here >>> + if (isinstance(object, (statesp.StateSpace, >>> xferfcn.TransferFunction))): >>> + return True >>> + >>> + # Didn't find anything that matched >>> + return False >>> >>> Modified: trunk/src/freqplot.py >>> =================================================================== >>> --- trunk/src/freqplot.py 2010-06-01 02:08:20 UTC (rev 20) >>> +++ trunk/src/freqplot.py 2010-06-06 17:29:42 UTC (rev 21) >>> @@ -42,11 +42,12 @@ >>> >>> import matplotlib.pyplot as plt >>> import scipy as sp >>> +import numpy as np >>> from ctrlutil import unwrap >>> from bdalg import feedback >>> >>> # Bode plot >>> -def bode(sys, omega=None, dB=False, Hz=False): >>> +def bode(syslist, omega=None, dB=False, Hz=False): >>> """Bode plot for a system >>> >>> Usage >>> @@ -57,8 +58,8 @@ >>> >>> Parameters >>> ---------- >>> - sys : linsys >>> - Linear input/output system >>> + syslist : linsys >>> + List of linear input/output systems (single system is OK) >>> omega : freq_range >>> Range of frequencies (list or bounds) in rad/sec >>> dB : boolean >>> @@ -76,45 +77,80 @@ >>> 1. Use (mag, phase, freq) = sys.freqresp(freq) to generate the >>> frequency response for a system. >>> """ >>> + # If argument was a singleton, turn it into a list >>> + if (not getattr(syslist, '__iter__', False)): >>> + syslist = (syslist,) >>> + >>> + # >>> # Select a default range if none is provided >>> + # >>> + # This code looks at the poles and zeros of all of the systems >>> that >>> + # we are plotting and sets the frequency range to be one >>> decade above >>> + # and below the min and max feature frequencies, rounded to >>> the nearest >>> + # integer. It excludes poles and zeros at the origin. If no >>> features >>> + # are found, it turns logspace(-1, 1) >>> + # >>> if (omega == None): >>> - omega = sp.logspace(-2, 2); >>> + # Find the list of all poles and zeros in the systems >>> + features = np.array(()) >>> + for sys in syslist: >>> + # Add new features to the list >>> + features = np.concatenate((features, np.abs(sys.poles))) >>> + features = np.concatenate((features, np.abs(sys.zeros))) >>> >>> - # Get the magnitude and phase of the system >>> - mag, phase, omega = sys.freqresp(omega) >>> - if Hz: >>> - omega = omega/(2*sp.pi) >>> - if dB: >>> - mag = 20*sp.log10(mag) >>> - phase = unwrap(phase*180/sp.pi, 360) >>> + # Get rid of poles and zeros at the origin >>> + features = features[features != 0]; >>> >>> - # Get the dimensions of the current axis, which we will divide >>> up >>> - #! TODO: Not current implemented; just use subplot for now >>> + # Make sure there is at least one point in the range >>> + if (features.shape[0] == 0): features = [1]; >>> >>> - # Magnitude plot >>> - plt.subplot(211); >>> - if dB: >>> - plt.semilogx(omega, mag) >>> - else: >>> - plt.loglog(omega, mag) >>> - plt.grid(True) >>> - plt.grid(True, which='minor') >>> - if dB: >>> - plt.ylabel("Magnitude (dB)") >>> - else: >>> - plt.ylabel("Magnitude") >>> + # Take the log of the features >>> + features = np.log10(features) >>> >>> - # Phase plot >>> - plt.subplot(212); >>> - plt.semilogx(omega, phase) >>> - plt.grid(True) >>> - plt.grid(True, which='minor') >>> - plt.ylabel("Phase (deg)") >>> - if Hz: >>> - plt.xlabel("Frequency (Hz)") >>> - else: >>> - plt.xlabel("Frequency (rad/sec)") >>> + # Set the to be an order of magnitude beyond any features >>> + omega = sp.logspace(np.floor(np.min(features))-1, >>> + np.ceil(np.max(features))+1) >>> >>> + for sys in syslist: >>> + # Get the magnitude and phase of the system >>> + mag, phase, omega = sys.freqresp(omega) >>> + if Hz: omega = omega/(2*sp.pi) >>> + if dB: mag = 20*sp.log10(mag) >>> + phase = unwrap(phase*180/sp.pi, 360) >>> + >>> + # Get the dimensions of the current axis, which we will >>> divide up >>> + #! TODO: Not current implemented; just use subplot for now >>> + >>> + # Magnitude plot >>> + plt.subplot(211); >>> + if dB: >>> + plt.semilogx(omega, mag) >>> + plt.ylabel("Magnitude (dB)") >>> + else: >>> + plt.loglog(omega, mag) >>> + plt.ylabel("Magnitude") >>> + >>> + # Add a grid to the plot >>> + plt.grid(True) >>> + plt.grid(True, which='minor') >>> + plt.hold(True); >>> + >>> + # Phase plot >>> + plt.subplot(212); >>> + plt.semilogx(omega, phase) >>> + plt.hold(True) >>> + >>> + # Add a grid to the plot >>> + plt.grid(True) >>> + plt.grid(True, which='minor') >>> + plt.ylabel("Phase (deg)") >>> + >>> + # Label the frequency axis >>> + if Hz: >>> + plt.xlabel("Frequency (Hz)") >>> + else: >>> + plt.xlabel("Frequency (rad/sec)") >>> + >>> return (211, 212) >>> >>> # Nyquist plot >>> >>> Modified: trunk/src/matlab.py >>> =================================================================== >>> --- trunk/src/matlab.py 2010-06-01 02:08:20 UTC (rev 20) >>> +++ trunk/src/matlab.py 2010-06-06 17:29:42 UTC (rev 21) >>> @@ -50,6 +50,7 @@ >>> >>> # Libraries that we make use of >>> import scipy as sp # SciPy library (used all over) >>> +import numpy as np # NumPy library >>> import scipy.signal as signal # Signal processing library >>> >>> # Import MATLAB-like functions that are defined in other packages >>> @@ -57,11 +58,16 @@ >>> from scipy.signal import lsim, impulse, step >>> from scipy import linspace, logspace >>> >>> -# Import MATLAB-like functions that belong elsewhere in python- >>> control >>> -from ctrlutil import unwrap >>> -from freqplot import bode, nyquist, gangof4 >>> +# Control system library >>> +import ctrlutil >>> +import freqplot >>> from statesp import StateSpace >>> from xferfcn import TransferFunction >>> +from exception import * >>> + >>> +# Import MATLAB-like functions that can be used as-is >>> +from ctrlutil import unwrap >>> +from freqplot import nyquist, gangof4 >>> from bdalg import series, parallel, negate, feedback >>> from pzmap import pzmap >>> from statefbk import place, lqr >>> @@ -304,3 +310,63 @@ >>> def freqresp(H, omega): >>> """Return the frequency response for an object H at frequency >>> omega""" >>> return H.freqresp(omega) >>> + >>> +# Bode plots >>> +def bode(*args, **keywords): >>> + """Bode plot of the frequency response >>> + >>> + Usage >>> + ===== >>> + bode(sys) >>> + bode(sys, w) >>> + bode(sys1, sys2, ..., sysN) >>> + bode(sys1, sys2, ..., sysN, w) >>> + bode(sys1, 'plotstyle1', ..., sysN, 'plotstyleN') >>> + """ >>> + >>> + # If the first argument is a list, then assume python-control >>> calling format >>> + if (getattr(args[0], '__iter__', False)): >>> + return freqplot.bode(*args, **keywords) >>> + >>> + # Otherwise, run through the arguments and collect up arguments >>> + syslist = []; plotstyle=[]; omega=None; >>> + i = 0; >>> + while i < len(args): >>> + # Check to see if this is a system of some sort >>> + if (ctrlutil.issys(args[i])): >>> + # Append the system to our list of systems >>> + syslist.append(args[i]) >>> + i += 1 >>> + >>> + # See if the next object is a plotsytle (string) >>> + if (i < len(args) and isinstance(args[i], str)): >>> + plotstyle.append(args[i]) >>> + i += 1 >>> + >>> + # Go on to the next argument >>> + continue >>> + >>> + # See if this is a frequency list >>> + elif (isinstance(args[i], (list, np.ndarray))): >>> + omega = args[i] >>> + i += 1 >>> + break >>> + >>> + else: >>> + raise ControlArgument("unrecognized argument type") >>> + >>> + # Check to make sure that we processed all arguments >>> + if (i < len(args)): >>> + raise ControlArgument("not all arguments processed") >>> + >>> + # Check to make sure we got the same number of plotstyles as >>> systems >>> + if (len(plotstyle) != 0 and len(syslist) != len(plotstyle)): >>> + raise ControlArgument("number of systems and plotstyles >>> should be equal") >>> + >>> + # Warn about unimplemented plotstyles >>> + #! TODO: remove this when plot styles are implemented in bode() >>> + if (len(plotstyle) != 0): >>> + print("Warning (matabl.bode): plot styles not implemented"); >>> + >>> + # Call the bode command >>> + return freqplot.bode(syslist, omega, **keywords) >>> >>> Modified: trunk/src/pzmap.py >>> =================================================================== >>> --- trunk/src/pzmap.py 2010-06-01 02:08:20 UTC (rev 20) >>> +++ trunk/src/pzmap.py 2010-06-06 17:29:42 UTC (rev 21) >>> @@ -45,6 +45,7 @@ >>> import xferfcn >>> >>> # Compute poles and zeros for a system >>> +#! TODO: extend this to handle state space systems >>> def pzmap(sys, Plot=True): >>> """Plot a pole/zero map for a transfer function""" >>> if (isinstance(sys, xferfcn.TransferFunction)): >>> >>> Modified: trunk/src/statesp.py >>> =================================================================== >>> --- trunk/src/statesp.py 2010-06-01 02:08:20 UTC (rev 20) >>> +++ trunk/src/statesp.py 2010-06-06 17:29:42 UTC (rev 21) >>> @@ -98,6 +98,18 @@ >>> #! TODO: Not implemented >>> return None >>> >>> + # Compute poles and zeros >>> + def poles(self): return sp.roots(sp.poly(self.A)) >>> + def zeros(self): >>> + den = sp.poly1d(sp.poly(self.A)) >>> + >>> + # Compute the numerator based on zeros >>> + #! TODO: This is currently limited to SISO systems >>> + num = sp.poly1d(\ >>> + sp.poly(self.A - sp.dot(self.B, self.C)) + (self.D[0] >>> - 1) * den) >>> + >>> + return (sp.roots(num)) >>> + >>> # Negation of a system >>> def __neg__(self): >>> """Negate a state space system""" >>> >>> Modified: trunk/src/xferfcn.py >>> =================================================================== >>> --- trunk/src/xferfcn.py 2010-06-01 02:08:20 UTC (rev 20) >>> +++ trunk/src/xferfcn.py 2010-06-06 17:29:42 UTC (rev 21) >>> @@ -190,6 +190,10 @@ >>> >>> return mag, phase, omega >>> >>> + # Compute poles and zeros >>> + def poles(self): return sp.roots(self.den) >>> + def zeros(self): return sp.roots(self.num) >>> + >>> # Feedback around a transfer function >>> def feedback(sys1, sys2, sign=-1): >>> """Feedback interconnection between two transfer functions""" >>> >>> >>> This was sent by the SourceForge.net collaborative development >>> platform, the world's largest Open Source development site. >> >> >> --- >> --- >> --- >> --------------------------------------------------------------------- >> ThinkGeek and WIRED's GeekDad team up for the Ultimate >> GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the >> lucky parental unit. See the prize list and enter to win: >> http://p.sf.net/sfu/thinkgeek-promo >> _______________________________________________ >> python-control-discuss mailing list >> pyt...@li... >> https://lists.sourceforge.net/lists/listinfo/python-control-discuss > > ------------------------------------------------------------------------------ > ThinkGeek and WIRED's GeekDad team up for the Ultimate > GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the > lucky parental unit. See the prize list and enter to win: > http://p.sf.net/sfu/thinkgeek-promo > _______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss |
|
From: <mur...@us...> - 2010-06-06 21:24:13
|
Revision: 22
http://python-control.svn.sourceforge.net/python-control/?rev=22&view=rev
Author: murrayrm
Date: 2010-06-06 21:24:07 +0000 (Sun, 06 Jun 2010)
Log Message:
-----------
preliminary LQR example (just MATLAB version for now)
Added Paths:
-----------
trunk/examples/pvtol-lqr.py
Added: trunk/examples/pvtol-lqr.py
===================================================================
--- trunk/examples/pvtol-lqr.py (rev 0)
+++ trunk/examples/pvtol-lqr.py 2010-06-06 21:24:07 UTC (rev 22)
@@ -0,0 +1,192 @@
+% pvtol_lqr.m - LQR design for vectored thrust aircraft
+% RMM, 14 Jan 03
+
+aminit;
+
+%%
+%% System dynamics
+%%
+%% These are the dynamics for the PVTOL system, written in state space
+%% form.
+%%
+
+pvtol_params; % System parameters
+
+% System matrices (entire plant: 2 input, 2 output)
+xe = [0 0 0 0 0 0]; ue = [0 m*g];
+[A, B, C, D] = pvtol_linearize(xe, ue);
+
+%%
+%% Construct inputs and outputs corresponding to steps in xy position
+%%
+%% The vectors xd and yd correspond to the states that are the desired
+%% equilibrium states for the system. The matrices Cx and Cy are the
+%% corresponding outputs.
+%%
+%% The way these vectors are used is to compute the closed loop system
+%% dynamics as
+%%
+%% xdot = Ax + B u => xdot = (A-BK)x + K xd
+%% u = -K(x - xd) y = Cx
+%%
+%% The closed loop dynamics can be simulated using the "step" command,
+%% with K*xd as the input vector (assumes that the "input" is unit size,
+%% so that xd corresponds to the desired steady state.
+%%
+
+xd = [1; 0; 0; 0; 0; 0]; Cx = [1 0 0 0 0 0];
+yd = [0; 1; 0; 0; 0; 0]; Cy = [0 1 0 0 0 0];
+
+%%
+%% LQR design
+%%
+
+% Start with a diagonal weighting
+Qx1 = diag([1, 1, 1, 1, 1, 1]);
+Qu1a = diag([1, 1]);
+K1a = lqr(A, B, Qx1, Qu1a);
+
+% Close the loop: xdot = Ax + B K (x-xd)
+H1a = ss(A-B*K1a, B*K1a*[xd, yd], [Cx; Cy], 0);
+[Y, T] = step(H1a, 10);
+
+figure(1); clf; subplot(321);
+ plot(T, Y(:,1, 1), '-', T, Y(:,2, 2), '--', ...
+ 'Linewidth', AM_data_linewidth); hold on;
+ plot([0 10], [1 1], 'k-', 'Linewidth', AM_ref_linewidth); hold on;
+
+ amaxis([0, 10, -0.1, 1.4]);
+ xlabel('time'); ylabel('position');
+
+ lgh = legend('x', 'y', 'Location', 'southeast');
+ legend(lgh, 'boxoff');
+ amprint('pvtol-lqrstep1.eps');
+
+% Look at different input weightings
+Qu1a = diag([1, 1]); K1a = lqr(A, B, Qx1, Qu1a);
+H1ax = ss(A-B*K1a,B(:,1)*K1a(1,:)*xd,Cx,0);
+
+Qu1b = 40^2*diag([1, 1]); K1b = lqr(A, B, Qx1, Qu1b);
+H1bx = ss(A-B*K1b,B(:,1)*K1b(1,:)*xd,Cx,0);
+
+Qu1c = 200^2*diag([1, 1]); K1c = lqr(A, B, Qx1, Qu1c);
+H1cx = ss(A-B*K1c,B(:,1)*K1c(1,:)*xd,Cx,0);
+
+[Y1, T1] = step(H1ax, 10);
+[Y2, T2] = step(H1bx, 10);
+[Y3, T3] = step(H1cx, 10);
+
+figure(2); clf; subplot(321);
+ plot(T1, Y1, 'b-', 'Linewidth', AM_data_linewidth); hold on;
+ plot(T2, Y2, 'b-', 'Linewidth', AM_data_linewidth); hold on;
+ plot(T3, Y3, 'b-', 'Linewidth', AM_data_linewidth); hold on;
+ plot([0 10], [1 1], 'k-', 'Linewidth', AM_ref_linewidth); hold on;
+
+ amaxis([0, 10, -0.1, 1.4]);
+ xlabel('time'); ylabel('position');
+
+ arcarrow([1.3 0.8], [5 0.45], -6);
+ text(5.3, 0.4, 'rho');
+
+ amprint('pvtol-lqrstep2.eps');
+
+% Output weighting - change Qx to use outputs
+Qx2 = [Cx; Cy]' * [Cx; Cy];
+Qu2 = 0.1 * diag([1, 1]);
+K2 = lqr(A, B, Qx2, Qu2);
+
+H2x = ss(A-B*K2,B(:,1)*K2(1,:)*xd,Cx,0);
+H2y = ss(A-B*K2,B(:,2)*K2(2,:)*yd,Cy,0);
+
+figure(3); step(H2x, H2y, 10);
+legend('x', 'y');
+
+%%
+%% Physically motivated weighting
+%%
+%% Shoot for 1 cm error in x, 10 cm error in y. Try to keep the angle
+%% less than 5 degrees in making the adjustments. Penalize side forces
+%% due to loss in efficiency.
+%%
+
+Qx3 = diag([100, 10, 2*pi/5, 0, 0, 0]);
+Qu3 = 0.1 * diag([1, 10]);
+K3 = lqr(A, B, Qx3, Qu3);
+
+H3x = ss(A-B*K3,B(:,1)*K3(1,:)*xd,Cx,0);
+H3y = ss(A-B*K3,B(:,2)*K3(2,:)*yd,Cy,0);
+figure(4); clf; subplot(221);
+step(H3x, H3y, 10);
+legend('x', 'y');
+
+%%
+%% Velocity control
+%%
+%% In this example, we modify the system so that we control the
+%% velocity of the system in the x direction. We ignore the
+%% dynamics in the vertical (y) direction. These dynamics demonstrate
+%% the role of the feedforward system since the equilibrium point
+%% corresponding to vd neq 0 requires a nonzero input.
+%%
+%% For this example, we use a control law u = -K(x-xd) + ud and convert
+%% this to the form u = -K x + N r, where r is the reference input and
+%% N is computed as described in class.
+%%
+
+% Extract system dynamics: theta, xdot, thdot
+Av = A([3 4 6], [3 4 6]);
+Bv = B([3 4 6], 1);
+Cv = [0 1 0]; % choose vx as output
+Dv = 0;
+
+% Design the feedback term using LQR
+Qxv = diag([2*pi/5, 10, 0]);
+Quv = 0.1;
+Kv = lqr(Av, Bv, Qxv, Quv);
+
+% Design the feedforward term by solve for eq pt in terms of reference r
+T = [Av Bv; Cv Dv]; % system matrix
+Nxu = T \ [0; 0; 0; 1]; % compute [Nx; Nu]
+Nx = Nxu(1:3); Nu = Nxu(4); % extract Nx and Nu
+N = Nu + Kv*Nx; % compute feedforward term
+
+%%
+%% Design #1: no feedforward input, ud
+%%
+
+Nv1 = [0; 1; 0];
+Hv1 = ss(Av-Bv*Kv, Bv*Kv*Nx, Cv, 0);
+step(Hv1, 10);
+
+%%
+%% Design #2: compute feedforward gain corresponding to equilibrium point
+%%
+
+Hv2 = ss(Av-Bv*Kv, Bv*N, Cv, 0);
+step(Hv2, 10);
+
+%%
+%% Design #3: integral action
+%%
+%% Add a new state to the system that is given by xidot = v - vd. We
+%% construct the control law by computing an LQR gain for the augmented
+%% system.
+%%
+
+Ai = [Av, [0; 0; 0]; [Cv, 0]];
+Bi = [Bv; 0];
+Ci = [Cv, 0];
+Di = Dv;
+
+% Design the feedback term, including weight on integrator error
+Qxi = diag([2*pi/5, 10, 0, 10]);
+Qui = 0.1;
+Ki = lqr(Ai, Bi, Qxi, Qui);
+
+% Desired state (augmented)
+xid = [0; 1; 0; 0];
+
+% Construct the closed loop system (including integrator)
+Hi = ss(Ai-Bi*Ki,Bi*Ki*xid - [0; 0; 0; Ci*xid],Ci,0);
+step(Hi, 10);
+
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: Tony V. <ton...@al...> - 2010-06-06 21:11:49
|
Doesn't look like I received it. I also checked my spam folder just in case. -tv On Jun 6, 2010, at 10:31, Richard Murray <mu...@cd...> wrote: > Tony and Erik: > > Can you let me know if you got the note below from SourceForge. I > just posted some changes to the package (needed them for an example > I'm putting together) and want to make sure that everyone is getting > copies of the changes that got posted. If not, I'll sort out what > we need to change to make sure you get these sorts of e-mails. > > Thanks. > > -richard > > Begin forwarded message: > >> From: mur...@us... >> Date: 6 June 2010 10:29:43 PDT >> To: mu...@cd... >> Subject: SF.net SVN: python-control:[21] trunk >> >> Revision: 21 >> http://python-control.svn.sourceforge.net/python-control/?rev=21&view=rev >> Author: murrayrm >> Date: 2010-06-06 17:29:42 +0000 (Sun, 06 Jun 2010) >> >> Log Message: >> ----------- >> Added changes to allow bode() to plot frequency resp for multiple >> systems >> * freqplot.bode now allows a list as the first argument >> * Added matlab.bode wrapper for MATLAB syntax: bode(sys1, ..., sysN) >> >> Added automatic frequency range determination in bode >> * Computes range of frequencies based on pole/zero locations >> * Eventually need to split this out as a ctrlutil function for use >> elsewhere >> >> Modified Paths: >> -------------- >> trunk/ChangeLog >> trunk/src/ctrlutil.py >> trunk/src/freqplot.py >> trunk/src/matlab.py >> trunk/src/pzmap.py >> trunk/src/statesp.py >> trunk/src/xferfcn.py >> >> Modified: trunk/ChangeLog >> =================================================================== >> --- trunk/ChangeLog 2010-06-01 02:08:20 UTC (rev 20) >> +++ trunk/ChangeLog 2010-06-06 17:29:42 UTC (rev 21) >> @@ -1,3 +1,25 @@ >> +2010-06-06 Richard Murray <murray@sumatra.local> >> + >> + * src/matlab.py (bode): created a wrapper that allows MATLAB >> style >> + arguments for bode (eg, bode(sys1, sys2)) >> + >> + * src/ctrlutil.py (issys): added new function to check if an >> object >> + is a system (state space or transfer function). Will generalize >> + this latter to look for other compatible classes. >> + >> + * src/freqplot.py (bode): Compute frequency range of bode plot >> based >> + on poles and zeros >> + (bode): Allow bode plot to be passed a list (or tuple) as the >> first >> + argument, in which case multiple bode plots are generated >> + >> + * src/statesp.py (StateSpace.zeros): new function to compute >> zeros >> + for a state space system >> + (StateSpace): defined new functions to compute poles of a state >> + space system >> + >> + * src/xferfcn.py (TransferFunction): defined new functions to >> + compute poles and zeros of a transfer function. >> + >> 2010-05-31 Richard Murray <murray@sumatra.local> >> >> * src/exception.py (ControlNotImplemented): added new exception, >> to >> >> Modified: trunk/src/ctrlutil.py >> =================================================================== >> --- trunk/src/ctrlutil.py 2010-06-01 02:08:20 UTC (rev 20) >> +++ trunk/src/ctrlutil.py 2010-06-06 17:29:42 UTC (rev 21) >> @@ -40,6 +40,12 @@ >> # >> # $Id$ >> >> +# Packages that we need access to >> +import scipy as sp >> +import statesp >> +import xferfcn >> + >> +# Specific functions that we use >> from scipy import pi >> >> # Utility function to unwrap an angle measurement >> @@ -78,3 +84,12 @@ >> >> # return the updated list >> return angle >> + >> +# Determine if an object is a system >> +def issys(object): >> + # Check for a member of one of the classes that we define here >> + if (isinstance(object, (statesp.StateSpace, >> xferfcn.TransferFunction))): >> + return True >> + >> + # Didn't find anything that matched >> + return False >> >> Modified: trunk/src/freqplot.py >> =================================================================== >> --- trunk/src/freqplot.py 2010-06-01 02:08:20 UTC (rev 20) >> +++ trunk/src/freqplot.py 2010-06-06 17:29:42 UTC (rev 21) >> @@ -42,11 +42,12 @@ >> >> import matplotlib.pyplot as plt >> import scipy as sp >> +import numpy as np >> from ctrlutil import unwrap >> from bdalg import feedback >> >> # Bode plot >> -def bode(sys, omega=None, dB=False, Hz=False): >> +def bode(syslist, omega=None, dB=False, Hz=False): >> """Bode plot for a system >> >> Usage >> @@ -57,8 +58,8 @@ >> >> Parameters >> ---------- >> - sys : linsys >> - Linear input/output system >> + syslist : linsys >> + List of linear input/output systems (single system is OK) >> omega : freq_range >> Range of frequencies (list or bounds) in rad/sec >> dB : boolean >> @@ -76,45 +77,80 @@ >> 1. Use (mag, phase, freq) = sys.freqresp(freq) to generate the >> frequency response for a system. >> """ >> + # If argument was a singleton, turn it into a list >> + if (not getattr(syslist, '__iter__', False)): >> + syslist = (syslist,) >> + >> + # >> # Select a default range if none is provided >> + # >> + # This code looks at the poles and zeros of all of the systems >> that >> + # we are plotting and sets the frequency range to be one >> decade above >> + # and below the min and max feature frequencies, rounded to >> the nearest >> + # integer. It excludes poles and zeros at the origin. If no >> features >> + # are found, it turns logspace(-1, 1) >> + # >> if (omega == None): >> - omega = sp.logspace(-2, 2); >> + # Find the list of all poles and zeros in the systems >> + features = np.array(()) >> + for sys in syslist: >> + # Add new features to the list >> + features = np.concatenate((features, np.abs(sys.poles))) >> + features = np.concatenate((features, np.abs(sys.zeros))) >> >> - # Get the magnitude and phase of the system >> - mag, phase, omega = sys.freqresp(omega) >> - if Hz: >> - omega = omega/(2*sp.pi) >> - if dB: >> - mag = 20*sp.log10(mag) >> - phase = unwrap(phase*180/sp.pi, 360) >> + # Get rid of poles and zeros at the origin >> + features = features[features != 0]; >> >> - # Get the dimensions of the current axis, which we will divide >> up >> - #! TODO: Not current implemented; just use subplot for now >> + # Make sure there is at least one point in the range >> + if (features.shape[0] == 0): features = [1]; >> >> - # Magnitude plot >> - plt.subplot(211); >> - if dB: >> - plt.semilogx(omega, mag) >> - else: >> - plt.loglog(omega, mag) >> - plt.grid(True) >> - plt.grid(True, which='minor') >> - if dB: >> - plt.ylabel("Magnitude (dB)") >> - else: >> - plt.ylabel("Magnitude") >> + # Take the log of the features >> + features = np.log10(features) >> >> - # Phase plot >> - plt.subplot(212); >> - plt.semilogx(omega, phase) >> - plt.grid(True) >> - plt.grid(True, which='minor') >> - plt.ylabel("Phase (deg)") >> - if Hz: >> - plt.xlabel("Frequency (Hz)") >> - else: >> - plt.xlabel("Frequency (rad/sec)") >> + # Set the to be an order of magnitude beyond any features >> + omega = sp.logspace(np.floor(np.min(features))-1, >> + np.ceil(np.max(features))+1) >> >> + for sys in syslist: >> + # Get the magnitude and phase of the system >> + mag, phase, omega = sys.freqresp(omega) >> + if Hz: omega = omega/(2*sp.pi) >> + if dB: mag = 20*sp.log10(mag) >> + phase = unwrap(phase*180/sp.pi, 360) >> + >> + # Get the dimensions of the current axis, which we will >> divide up >> + #! TODO: Not current implemented; just use subplot for now >> + >> + # Magnitude plot >> + plt.subplot(211); >> + if dB: >> + plt.semilogx(omega, mag) >> + plt.ylabel("Magnitude (dB)") >> + else: >> + plt.loglog(omega, mag) >> + plt.ylabel("Magnitude") >> + >> + # Add a grid to the plot >> + plt.grid(True) >> + plt.grid(True, which='minor') >> + plt.hold(True); >> + >> + # Phase plot >> + plt.subplot(212); >> + plt.semilogx(omega, phase) >> + plt.hold(True) >> + >> + # Add a grid to the plot >> + plt.grid(True) >> + plt.grid(True, which='minor') >> + plt.ylabel("Phase (deg)") >> + >> + # Label the frequency axis >> + if Hz: >> + plt.xlabel("Frequency (Hz)") >> + else: >> + plt.xlabel("Frequency (rad/sec)") >> + >> return (211, 212) >> >> # Nyquist plot >> >> Modified: trunk/src/matlab.py >> =================================================================== >> --- trunk/src/matlab.py 2010-06-01 02:08:20 UTC (rev 20) >> +++ trunk/src/matlab.py 2010-06-06 17:29:42 UTC (rev 21) >> @@ -50,6 +50,7 @@ >> >> # Libraries that we make use of >> import scipy as sp # SciPy library (used all over) >> +import numpy as np # NumPy library >> import scipy.signal as signal # Signal processing library >> >> # Import MATLAB-like functions that are defined in other packages >> @@ -57,11 +58,16 @@ >> from scipy.signal import lsim, impulse, step >> from scipy import linspace, logspace >> >> -# Import MATLAB-like functions that belong elsewhere in python- >> control >> -from ctrlutil import unwrap >> -from freqplot import bode, nyquist, gangof4 >> +# Control system library >> +import ctrlutil >> +import freqplot >> from statesp import StateSpace >> from xferfcn import TransferFunction >> +from exception import * >> + >> +# Import MATLAB-like functions that can be used as-is >> +from ctrlutil import unwrap >> +from freqplot import nyquist, gangof4 >> from bdalg import series, parallel, negate, feedback >> from pzmap import pzmap >> from statefbk import place, lqr >> @@ -304,3 +310,63 @@ >> def freqresp(H, omega): >> """Return the frequency response for an object H at frequency >> omega""" >> return H.freqresp(omega) >> + >> +# Bode plots >> +def bode(*args, **keywords): >> + """Bode plot of the frequency response >> + >> + Usage >> + ===== >> + bode(sys) >> + bode(sys, w) >> + bode(sys1, sys2, ..., sysN) >> + bode(sys1, sys2, ..., sysN, w) >> + bode(sys1, 'plotstyle1', ..., sysN, 'plotstyleN') >> + """ >> + >> + # If the first argument is a list, then assume python-control >> calling format >> + if (getattr(args[0], '__iter__', False)): >> + return freqplot.bode(*args, **keywords) >> + >> + # Otherwise, run through the arguments and collect up arguments >> + syslist = []; plotstyle=[]; omega=None; >> + i = 0; >> + while i < len(args): >> + # Check to see if this is a system of some sort >> + if (ctrlutil.issys(args[i])): >> + # Append the system to our list of systems >> + syslist.append(args[i]) >> + i += 1 >> + >> + # See if the next object is a plotsytle (string) >> + if (i < len(args) and isinstance(args[i], str)): >> + plotstyle.append(args[i]) >> + i += 1 >> + >> + # Go on to the next argument >> + continue >> + >> + # See if this is a frequency list >> + elif (isinstance(args[i], (list, np.ndarray))): >> + omega = args[i] >> + i += 1 >> + break >> + >> + else: >> + raise ControlArgument("unrecognized argument type") >> + >> + # Check to make sure that we processed all arguments >> + if (i < len(args)): >> + raise ControlArgument("not all arguments processed") >> + >> + # Check to make sure we got the same number of plotstyles as >> systems >> + if (len(plotstyle) != 0 and len(syslist) != len(plotstyle)): >> + raise ControlArgument("number of systems and plotstyles >> should be equal") >> + >> + # Warn about unimplemented plotstyles >> + #! TODO: remove this when plot styles are implemented in bode() >> + if (len(plotstyle) != 0): >> + print("Warning (matabl.bode): plot styles not implemented"); >> + >> + # Call the bode command >> + return freqplot.bode(syslist, omega, **keywords) >> >> Modified: trunk/src/pzmap.py >> =================================================================== >> --- trunk/src/pzmap.py 2010-06-01 02:08:20 UTC (rev 20) >> +++ trunk/src/pzmap.py 2010-06-06 17:29:42 UTC (rev 21) >> @@ -45,6 +45,7 @@ >> import xferfcn >> >> # Compute poles and zeros for a system >> +#! TODO: extend this to handle state space systems >> def pzmap(sys, Plot=True): >> """Plot a pole/zero map for a transfer function""" >> if (isinstance(sys, xferfcn.TransferFunction)): >> >> Modified: trunk/src/statesp.py >> =================================================================== >> --- trunk/src/statesp.py 2010-06-01 02:08:20 UTC (rev 20) >> +++ trunk/src/statesp.py 2010-06-06 17:29:42 UTC (rev 21) >> @@ -98,6 +98,18 @@ >> #! TODO: Not implemented >> return None >> >> + # Compute poles and zeros >> + def poles(self): return sp.roots(sp.poly(self.A)) >> + def zeros(self): >> + den = sp.poly1d(sp.poly(self.A)) >> + >> + # Compute the numerator based on zeros >> + #! TODO: This is currently limited to SISO systems >> + num = sp.poly1d(\ >> + sp.poly(self.A - sp.dot(self.B, self.C)) + (self.D[0] >> - 1) * den) >> + >> + return (sp.roots(num)) >> + >> # Negation of a system >> def __neg__(self): >> """Negate a state space system""" >> >> Modified: trunk/src/xferfcn.py >> =================================================================== >> --- trunk/src/xferfcn.py 2010-06-01 02:08:20 UTC (rev 20) >> +++ trunk/src/xferfcn.py 2010-06-06 17:29:42 UTC (rev 21) >> @@ -190,6 +190,10 @@ >> >> return mag, phase, omega >> >> + # Compute poles and zeros >> + def poles(self): return sp.roots(self.den) >> + def zeros(self): return sp.roots(self.num) >> + >> # Feedback around a transfer function >> def feedback(sys1, sys2, sign=-1): >> """Feedback interconnection between two transfer functions""" >> >> >> This was sent by the SourceForge.net collaborative development >> platform, the world's largest Open Source development site. > > > --- > --- > --- > --------------------------------------------------------------------- > ThinkGeek and WIRED's GeekDad team up for the Ultimate > GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the > lucky parental unit. See the prize list and enter to win: > http://p.sf.net/sfu/thinkgeek-promo > _______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss |
|
From: Richard M. <mu...@cd...> - 2010-06-06 17:31:32
|
Tony and Erik: Can you let me know if you got the note below from SourceForge. I just posted some changes to the package (needed them for an example I'm putting together) and want to make sure that everyone is getting copies of the changes that got posted. If not, I'll sort out what we need to change to make sure you get these sorts of e-mails. Thanks. -richard Begin forwarded message: > From: mur...@us... > Date: 6 June 2010 10:29:43 PDT > To: mu...@cd... > Subject: SF.net SVN: python-control:[21] trunk > > Revision: 21 > http://python-control.svn.sourceforge.net/python-control/?rev=21&view=rev > Author: murrayrm > Date: 2010-06-06 17:29:42 +0000 (Sun, 06 Jun 2010) > > Log Message: > ----------- > Added changes to allow bode() to plot frequency resp for multiple systems > * freqplot.bode now allows a list as the first argument > * Added matlab.bode wrapper for MATLAB syntax: bode(sys1, ..., sysN) > > Added automatic frequency range determination in bode > * Computes range of frequencies based on pole/zero locations > * Eventually need to split this out as a ctrlutil function for use elsewhere > > Modified Paths: > -------------- > trunk/ChangeLog > trunk/src/ctrlutil.py > trunk/src/freqplot.py > trunk/src/matlab.py > trunk/src/pzmap.py > trunk/src/statesp.py > trunk/src/xferfcn.py > > Modified: trunk/ChangeLog > =================================================================== > --- trunk/ChangeLog 2010-06-01 02:08:20 UTC (rev 20) > +++ trunk/ChangeLog 2010-06-06 17:29:42 UTC (rev 21) > @@ -1,3 +1,25 @@ > +2010-06-06 Richard Murray <murray@sumatra.local> > + > + * src/matlab.py (bode): created a wrapper that allows MATLAB style > + arguments for bode (eg, bode(sys1, sys2)) > + > + * src/ctrlutil.py (issys): added new function to check if an object > + is a system (state space or transfer function). Will generalize > + this latter to look for other compatible classes. > + > + * src/freqplot.py (bode): Compute frequency range of bode plot based > + on poles and zeros > + (bode): Allow bode plot to be passed a list (or tuple) as the first > + argument, in which case multiple bode plots are generated > + > + * src/statesp.py (StateSpace.zeros): new function to compute zeros > + for a state space system > + (StateSpace): defined new functions to compute poles of a state > + space system > + > + * src/xferfcn.py (TransferFunction): defined new functions to > + compute poles and zeros of a transfer function. > + > 2010-05-31 Richard Murray <murray@sumatra.local> > > * src/exception.py (ControlNotImplemented): added new exception, to > > Modified: trunk/src/ctrlutil.py > =================================================================== > --- trunk/src/ctrlutil.py 2010-06-01 02:08:20 UTC (rev 20) > +++ trunk/src/ctrlutil.py 2010-06-06 17:29:42 UTC (rev 21) > @@ -40,6 +40,12 @@ > # > # $Id$ > > +# Packages that we need access to > +import scipy as sp > +import statesp > +import xferfcn > + > +# Specific functions that we use > from scipy import pi > > # Utility function to unwrap an angle measurement > @@ -78,3 +84,12 @@ > > # return the updated list > return angle > + > +# Determine if an object is a system > +def issys(object): > + # Check for a member of one of the classes that we define here > + if (isinstance(object, (statesp.StateSpace, xferfcn.TransferFunction))): > + return True > + > + # Didn't find anything that matched > + return False > > Modified: trunk/src/freqplot.py > =================================================================== > --- trunk/src/freqplot.py 2010-06-01 02:08:20 UTC (rev 20) > +++ trunk/src/freqplot.py 2010-06-06 17:29:42 UTC (rev 21) > @@ -42,11 +42,12 @@ > > import matplotlib.pyplot as plt > import scipy as sp > +import numpy as np > from ctrlutil import unwrap > from bdalg import feedback > > # Bode plot > -def bode(sys, omega=None, dB=False, Hz=False): > +def bode(syslist, omega=None, dB=False, Hz=False): > """Bode plot for a system > > Usage > @@ -57,8 +58,8 @@ > > Parameters > ---------- > - sys : linsys > - Linear input/output system > + syslist : linsys > + List of linear input/output systems (single system is OK) > omega : freq_range > Range of frequencies (list or bounds) in rad/sec > dB : boolean > @@ -76,45 +77,80 @@ > 1. Use (mag, phase, freq) = sys.freqresp(freq) to generate the > frequency response for a system. > """ > + # If argument was a singleton, turn it into a list > + if (not getattr(syslist, '__iter__', False)): > + syslist = (syslist,) > + > + # > # Select a default range if none is provided > + # > + # This code looks at the poles and zeros of all of the systems that > + # we are plotting and sets the frequency range to be one decade above > + # and below the min and max feature frequencies, rounded to the nearest > + # integer. It excludes poles and zeros at the origin. If no features > + # are found, it turns logspace(-1, 1) > + # > if (omega == None): > - omega = sp.logspace(-2, 2); > + # Find the list of all poles and zeros in the systems > + features = np.array(()) > + for sys in syslist: > + # Add new features to the list > + features = np.concatenate((features, np.abs(sys.poles))) > + features = np.concatenate((features, np.abs(sys.zeros))) > > - # Get the magnitude and phase of the system > - mag, phase, omega = sys.freqresp(omega) > - if Hz: > - omega = omega/(2*sp.pi) > - if dB: > - mag = 20*sp.log10(mag) > - phase = unwrap(phase*180/sp.pi, 360) > + # Get rid of poles and zeros at the origin > + features = features[features != 0]; > > - # Get the dimensions of the current axis, which we will divide up > - #! TODO: Not current implemented; just use subplot for now > + # Make sure there is at least one point in the range > + if (features.shape[0] == 0): features = [1]; > > - # Magnitude plot > - plt.subplot(211); > - if dB: > - plt.semilogx(omega, mag) > - else: > - plt.loglog(omega, mag) > - plt.grid(True) > - plt.grid(True, which='minor') > - if dB: > - plt.ylabel("Magnitude (dB)") > - else: > - plt.ylabel("Magnitude") > + # Take the log of the features > + features = np.log10(features) > > - # Phase plot > - plt.subplot(212); > - plt.semilogx(omega, phase) > - plt.grid(True) > - plt.grid(True, which='minor') > - plt.ylabel("Phase (deg)") > - if Hz: > - plt.xlabel("Frequency (Hz)") > - else: > - plt.xlabel("Frequency (rad/sec)") > + # Set the to be an order of magnitude beyond any features > + omega = sp.logspace(np.floor(np.min(features))-1, > + np.ceil(np.max(features))+1) > > + for sys in syslist: > + # Get the magnitude and phase of the system > + mag, phase, omega = sys.freqresp(omega) > + if Hz: omega = omega/(2*sp.pi) > + if dB: mag = 20*sp.log10(mag) > + phase = unwrap(phase*180/sp.pi, 360) > + > + # Get the dimensions of the current axis, which we will divide up > + #! TODO: Not current implemented; just use subplot for now > + > + # Magnitude plot > + plt.subplot(211); > + if dB: > + plt.semilogx(omega, mag) > + plt.ylabel("Magnitude (dB)") > + else: > + plt.loglog(omega, mag) > + plt.ylabel("Magnitude") > + > + # Add a grid to the plot > + plt.grid(True) > + plt.grid(True, which='minor') > + plt.hold(True); > + > + # Phase plot > + plt.subplot(212); > + plt.semilogx(omega, phase) > + plt.hold(True) > + > + # Add a grid to the plot > + plt.grid(True) > + plt.grid(True, which='minor') > + plt.ylabel("Phase (deg)") > + > + # Label the frequency axis > + if Hz: > + plt.xlabel("Frequency (Hz)") > + else: > + plt.xlabel("Frequency (rad/sec)") > + > return (211, 212) > > # Nyquist plot > > Modified: trunk/src/matlab.py > =================================================================== > --- trunk/src/matlab.py 2010-06-01 02:08:20 UTC (rev 20) > +++ trunk/src/matlab.py 2010-06-06 17:29:42 UTC (rev 21) > @@ -50,6 +50,7 @@ > > # Libraries that we make use of > import scipy as sp # SciPy library (used all over) > +import numpy as np # NumPy library > import scipy.signal as signal # Signal processing library > > # Import MATLAB-like functions that are defined in other packages > @@ -57,11 +58,16 @@ > from scipy.signal import lsim, impulse, step > from scipy import linspace, logspace > > -# Import MATLAB-like functions that belong elsewhere in python-control > -from ctrlutil import unwrap > -from freqplot import bode, nyquist, gangof4 > +# Control system library > +import ctrlutil > +import freqplot > from statesp import StateSpace > from xferfcn import TransferFunction > +from exception import * > + > +# Import MATLAB-like functions that can be used as-is > +from ctrlutil import unwrap > +from freqplot import nyquist, gangof4 > from bdalg import series, parallel, negate, feedback > from pzmap import pzmap > from statefbk import place, lqr > @@ -304,3 +310,63 @@ > def freqresp(H, omega): > """Return the frequency response for an object H at frequency omega""" > return H.freqresp(omega) > + > +# Bode plots > +def bode(*args, **keywords): > + """Bode plot of the frequency response > + > + Usage > + ===== > + bode(sys) > + bode(sys, w) > + bode(sys1, sys2, ..., sysN) > + bode(sys1, sys2, ..., sysN, w) > + bode(sys1, 'plotstyle1', ..., sysN, 'plotstyleN') > + """ > + > + # If the first argument is a list, then assume python-control calling format > + if (getattr(args[0], '__iter__', False)): > + return freqplot.bode(*args, **keywords) > + > + # Otherwise, run through the arguments and collect up arguments > + syslist = []; plotstyle=[]; omega=None; > + i = 0; > + while i < len(args): > + # Check to see if this is a system of some sort > + if (ctrlutil.issys(args[i])): > + # Append the system to our list of systems > + syslist.append(args[i]) > + i += 1 > + > + # See if the next object is a plotsytle (string) > + if (i < len(args) and isinstance(args[i], str)): > + plotstyle.append(args[i]) > + i += 1 > + > + # Go on to the next argument > + continue > + > + # See if this is a frequency list > + elif (isinstance(args[i], (list, np.ndarray))): > + omega = args[i] > + i += 1 > + break > + > + else: > + raise ControlArgument("unrecognized argument type") > + > + # Check to make sure that we processed all arguments > + if (i < len(args)): > + raise ControlArgument("not all arguments processed") > + > + # Check to make sure we got the same number of plotstyles as systems > + if (len(plotstyle) != 0 and len(syslist) != len(plotstyle)): > + raise ControlArgument("number of systems and plotstyles should be equal") > + > + # Warn about unimplemented plotstyles > + #! TODO: remove this when plot styles are implemented in bode() > + if (len(plotstyle) != 0): > + print("Warning (matabl.bode): plot styles not implemented"); > + > + # Call the bode command > + return freqplot.bode(syslist, omega, **keywords) > > Modified: trunk/src/pzmap.py > =================================================================== > --- trunk/src/pzmap.py 2010-06-01 02:08:20 UTC (rev 20) > +++ trunk/src/pzmap.py 2010-06-06 17:29:42 UTC (rev 21) > @@ -45,6 +45,7 @@ > import xferfcn > > # Compute poles and zeros for a system > +#! TODO: extend this to handle state space systems > def pzmap(sys, Plot=True): > """Plot a pole/zero map for a transfer function""" > if (isinstance(sys, xferfcn.TransferFunction)): > > Modified: trunk/src/statesp.py > =================================================================== > --- trunk/src/statesp.py 2010-06-01 02:08:20 UTC (rev 20) > +++ trunk/src/statesp.py 2010-06-06 17:29:42 UTC (rev 21) > @@ -98,6 +98,18 @@ > #! TODO: Not implemented > return None > > + # Compute poles and zeros > + def poles(self): return sp.roots(sp.poly(self.A)) > + def zeros(self): > + den = sp.poly1d(sp.poly(self.A)) > + > + # Compute the numerator based on zeros > + #! TODO: This is currently limited to SISO systems > + num = sp.poly1d(\ > + sp.poly(self.A - sp.dot(self.B, self.C)) + (self.D[0] - 1) * den) > + > + return (sp.roots(num)) > + > # Negation of a system > def __neg__(self): > """Negate a state space system""" > > Modified: trunk/src/xferfcn.py > =================================================================== > --- trunk/src/xferfcn.py 2010-06-01 02:08:20 UTC (rev 20) > +++ trunk/src/xferfcn.py 2010-06-06 17:29:42 UTC (rev 21) > @@ -190,6 +190,10 @@ > > return mag, phase, omega > > + # Compute poles and zeros > + def poles(self): return sp.roots(self.den) > + def zeros(self): return sp.roots(self.num) > + > # Feedback around a transfer function > def feedback(sys1, sys2, sign=-1): > """Feedback interconnection between two transfer functions""" > > > This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
|
From: Tony V. <ton...@al...> - 2010-06-05 01:46:05
|
I agree with running against multiple versions; I do the same thing with some of the tools I build at JPL. So that's pretty much what I thought you'd say. Just thought I'd check. -tv On Jun 4, 2010, at 17:59, Richard Murray <mu...@cd...> wrote: > It would be good if the package itself ran on as many versions of > pythons as we can easily support. The CDS cluster is running python > 2.4 (about 2 years out of date). I've got python 2.5.4 and 2.6.1 > (via OS X) as well as python 2.6.5 (via fink) on my Mac. By > default, I've been running 2.6.1 for all of my testing. > > There is some advantage to running different version of python since > that will help us keep compatibility across multiple flavors. For > what we are doing, I don't think there should be anything version > specific that comes up. As you've probably found, it's mainly just > a matter of getting a version of ipython up an running. This is not > quite as easy as it should be, unfortunately... > > FWIW, I used fink for my ipython and scipy installs. Even then, I > had to turn on the unstable distribution. > > -richard > > On 4 Jun 2010, at 17:35 , Tony Vanelli wrote: > >> Richard, >> >> After a fair bit of wrestling with the fact that various >> ramifications of OS X *also* having a older version of Python than >> what I had downloaded and installed, I've gotten the the Python >> Control Library v.0.3b to execute the 2nd order test-case. (Not as >> far as I wanted to get, but it's a start...) >> >> But my experience of dealing with different Python versions >> (described below) on OS X has raised a question: >> >> Do you have a preferred version of Python you'd like us to use? >> >> OS X 10.6 comes bundled Python 2.6.1; I've installed Python 2.6.5 >> (latest non-developmental release). I recall Erik mentioned he'd >> installed Py2.7 in an earlier email. >> >> Erik: If you've installed Py2.7 and not tried running the 2nd-order >> example case yet, then you might be about to run into the same >> issues I smacked into. It turns out that various "default" >> packages (and tools) that either come with Python 2.6.1, or that >> Apple pre-builds for you, are *not* included in the Python 2.6.4+ >> distros available from python.org, and some these packages have to >> be installed by hand. >> >> iPython, it turns out, needs a little extra care. My first >> installation apparently "locked on" to the Apple Python distro >> rather than the one I installed. A simple first check is to run >> iPython (no arguments) and see which version of Python you *really* >> get... see excerpt below; note different boldfaced version numbers. >> >> I've got things cobbled together now, but I'm going to spend a >> little time ensuring it's clean. If you hit a snag, I (fairly >> soon) should have some tips I can offer. >> >> -tv >> >> 43 acheron: /Users/tonyv> ipython >> Python 2.6.1 (r261:67515, Feb 11 2010, 00:51:29) >> Type "copyright", "credits" or "license" for more information. >> >> IPython 0.10 -- An enhanced Interactive Python. >> ? -> Introduction and overview of IPython's features. >> %quickref -> Quick reference. >> help -> Python's own help system. >> object? -> Details about 'object'. ?object also works, ?? prints >> more. >> >> In [1]: quit() >> Do you really want to exit ([y]/n)? y >> >> 44 acheron: /Users/tonyv> python >> Python 2.6.5 (r265:79359, Mar 24 2010, 01:32:55) >> [GCC 4.0.1 (Apple Inc. build 5493)] on darwin >> Type "help", "copyright", "credits" or "license" for more >> information. >>>>> quit() >> 45 acheron: /Users/tonyv> >> >> >> >> --- >> --- >> --- >> --------------------------------------------------------------------- >> ThinkGeek and WIRED's GeekDad team up for the Ultimate >> GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the >> lucky parental unit. See the prize list and enter to win: >> http://p.sf.net/sfu/thinkgeek-promo_______________________________________________ >> python-control-discuss mailing list >> pyt...@li... >> https://lists.sourceforge.net/lists/listinfo/python-control-discuss > > > --- > --- > --- > --------------------------------------------------------------------- > ThinkGeek and WIRED's GeekDad team up for the Ultimate > GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the > lucky parental unit. See the prize list and enter to win: > http://p.sf.net/sfu/thinkgeek-promo > _______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss |
|
From: Richard M. <mu...@cd...> - 2010-06-05 00:59:20
|
It would be good if the package itself ran on as many versions of pythons as we can easily support. The CDS cluster is running python 2.4 (about 2 years out of date). I've got python 2.5.4 and 2.6.1 (via OS X) as well as python 2.6.5 (via fink) on my Mac. By default, I've been running 2.6.1 for all of my testing. There is some advantage to running different version of python since that will help us keep compatibility across multiple flavors. For what we are doing, I don't think there should be anything version specific that comes up. As you've probably found, it's mainly just a matter of getting a version of ipython up an running. This is not quite as easy as it should be, unfortunately... FWIW, I used fink for my ipython and scipy installs. Even then, I had to turn on the unstable distribution. -richard On 4 Jun 2010, at 17:35 , Tony Vanelli wrote: > Richard, > > After a fair bit of wrestling with the fact that various ramifications of OS X *also* having a older version of Python than what I had downloaded and installed, I've gotten the the Python Control Library v.0.3b to execute the 2nd order test-case. (Not as far as I wanted to get, but it's a start...) > > But my experience of dealing with different Python versions (described below) on OS X has raised a question: > > Do you have a preferred version of Python you'd like us to use? > > OS X 10.6 comes bundled Python 2.6.1; I've installed Python 2.6.5 (latest non-developmental release). I recall Erik mentioned he'd installed Py2.7 in an earlier email. > > Erik: If you've installed Py2.7 and not tried running the 2nd-order example case yet, then you might be about to run into the same issues I smacked into. It turns out that various "default" packages (and tools) that either come with Python 2.6.1, or that Apple pre-builds for you, are *not* included in the Python 2.6.4+ distros available from python.org, and some these packages have to be installed by hand. > > iPython, it turns out, needs a little extra care. My first installation apparently "locked on" to the Apple Python distro rather than the one I installed. A simple first check is to run iPython (no arguments) and see which version of Python you *really* get... see excerpt below; note different boldfaced version numbers. > > I've got things cobbled together now, but I'm going to spend a little time ensuring it's clean. If you hit a snag, I (fairly soon) should have some tips I can offer. > > -tv > > 43 acheron: /Users/tonyv> ipython > Python 2.6.1 (r261:67515, Feb 11 2010, 00:51:29) > Type "copyright", "credits" or "license" for more information. > > IPython 0.10 -- An enhanced Interactive Python. > ? -> Introduction and overview of IPython's features. > %quickref -> Quick reference. > help -> Python's own help system. > object? -> Details about 'object'. ?object also works, ?? prints more. > > In [1]: quit() > Do you really want to exit ([y]/n)? y > > 44 acheron: /Users/tonyv> python > Python 2.6.5 (r265:79359, Mar 24 2010, 01:32:55) > [GCC 4.0.1 (Apple Inc. build 5493)] on darwin > Type "help", "copyright", "credits" or "license" for more information. > >>> quit() > 45 acheron: /Users/tonyv> > > > > ------------------------------------------------------------------------------ > ThinkGeek and WIRED's GeekDad team up for the Ultimate > GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the > lucky parental unit. See the prize list and enter to win: > http://p.sf.net/sfu/thinkgeek-promo_______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss |
|
From: Tony V. <ton...@al...> - 2010-06-05 00:36:09
|
Richard, After a fair bit of wrestling with the fact that various ramifications of OS X *also* having a older version of Python than what I had downloaded and installed, I've gotten the the Python Control Library v.0.3b to execute the 2nd order test-case. (Not as far as I wanted to get, but it's a start...) But my experience of dealing with different Python versions (described below) on OS X has raised a question: Do you have a preferred version of Python you'd like us to use? OS X 10.6 comes bundled Python 2.6.1; I've installed Python 2.6.5 (latest non-developmental release). I recall Erik mentioned he'd installed Py2.7 in an earlier email. Erik: If you've installed Py2.7 and not tried running the 2nd-order example case yet, then you might be about to run into the same issues I smacked into. It turns out that various "default" packages (and tools) that either come with Python 2.6.1, or that Apple pre-builds for you, are *not* included in the Python 2.6.4+ distros available from python.org, and some these packages have to be installed by hand. iPython, it turns out, needs a little extra care. My first installation apparently "locked on" to the Apple Python distro rather than the one I installed. A simple first check is to run iPython (no arguments) and see which version of Python you *really* get... see excerpt below; note different boldfaced version numbers. I've got things cobbled together now, but I'm going to spend a little time ensuring it's clean. If you hit a snag, I (fairly soon) should have some tips I can offer. -tv 43 acheron: /Users/tonyv> ipython Python 2.6.1 (r261:67515, Feb 11 2010, 00:51:29) Type "copyright", "credits" or "license" for more information. IPython 0.10 -- An enhanced Interactive Python. ? -> Introduction and overview of IPython's features. %quickref -> Quick reference. help -> Python's own help system. object? -> Details about 'object'. ?object also works, ?? prints more. In [1]: quit() Do you really want to exit ([y]/n)? y 44 acheron: /Users/tonyv> python Python 2.6.5 (r265:79359, Mar 24 2010, 01:32:55) [GCC 4.0.1 (Apple Inc. build 5493)] on darwin Type "help", "copyright", "credits" or "license" for more information. >>> quit() 45 acheron: /Users/tonyv> |
|
From: Tony V. <ton...@al...> - 2010-06-03 06:35:10
|
Got'em. I'll try to review the state-of-things-so-far this Friday, per the earlier email notes you sent Erik. -tv On Jun 2, 2010, at 3:08 PM, Richard Murray wrote: > Erik and Tony, > > You should have just been subscribed to two separate lists: > > * python-control-discuss (this list) > * python-control-announce - announcements about the package > > The plan is to use this list to have back and forth discussions related to the package, so that everyone can see the discussion. If it gets to overwhelming (not likely, given that we are all doing this in our spare time!) we can use some of the other mechanisms on sourceforge (forums, tracking) to organize the discussions better. > > Kristian, Gunnar and Enrico: I signed you all up for the announce list since you have been interacting with me on the project in some way and I figured you might be interested in occasional e-mail as functionality gets added. If you want to be more involved in either development or discussions, I'm happy to add you to this list as well. Just let me know. > > -richard > > > ------------------------------------------------------------------------------ > ThinkGeek and WIRED's GeekDad team up for the Ultimate > GeekDad Father's Day Giveaway. ONE MASSIVE PRIZE to the > lucky parental unit. See the prize list and enter to win: > http://p.sf.net/sfu/thinkgeek-promo > _______________________________________________ > python-control-discuss mailing list > pyt...@li... > https://lists.sourceforge.net/lists/listinfo/python-control-discuss |
|
From: Richard M. <mu...@cd...> - 2010-06-02 22:08:23
|
Erik and Tony, You should have just been subscribed to two separate lists: * python-control-discuss (this list) * python-control-announce - announcements about the package The plan is to use this list to have back and forth discussions related to the package, so that everyone can see the discussion. If it gets to overwhelming (not likely, given that we are all doing this in our spare time!) we can use some of the other mechanisms on sourceforge (forums, tracking) to organize the discussions better. Kristian, Gunnar and Enrico: I signed you all up for the announce list since you have been interacting with me on the project in some way and I figured you might be interested in occasional e-mail as functionality gets added. If you want to be more involved in either development or discussions, I'm happy to add you to this list as well. Just let me know. -richard |
7 messages has been excluded from this view by a project administrator.
