Movatterモバイル変換

[0]ホーム
URL:

Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly

 Sign up
Appearance settings

Commitb72da8a

Browse files
committed
Added functionality to balred() in modelsimp.py. Now does stable and unstable systems along with 'matchdc' gain option. Also added option to gram() to return cholesky factored gramian.
1 parentb10a823 commitb72da8a

File tree

2 files changed

+105
-54
lines changed

2 files changed

+105
-54
lines changed

‎control/modelsimp.py

Lines changed: 67 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -197,6 +197,37 @@ def modred(sys, ELIM, method='matchdc'):
197197
rsys=StateSpace(Ar,Br,Cr,Dr)
198198
returnrsys
199199

200+
defstabsep(T_schur,Z_schur,sys,ldim,no_in,no_out):
201+
"""
202+
Performs stable/unstabe decomposition of sys after Schur forms have been computed for system matrix.
203+
204+
Reference: Hsu,C.S., and Hou,D., 1991,
205+
Reducing unstable linear control systems via real Schur transformation.
206+
Electronics Letters, 27, 984-986.
207+
208+
"""
209+
#Author: M. Clement (mdclemen@eng.ucsd.edu) 2016
210+
As=np.asmatrix(T_schur)
211+
Bs=Z_schur.T*sys.B
212+
Cs=sys.C*Z_schur
213+
#from ref 1 eq(1) As = [A_ Ac], Bs = [B_], and Cs = [C_ C+]; _ denotes stable subsystem
214+
# [0 A+] [B+]
215+
A_=As[0:ldim,0:ldim]
216+
Ac=As[0:ldim,ldim::]
217+
Ap=As[ldim::,ldim::]
218+
219+
B_=Bs[0:ldim,:]
220+
Bp=Bs[ldim::,:]
221+
222+
C_=Cs[:,0:ldim]
223+
Cp=Cs[:,ldim::]
224+
#do some more tricky math IAW ref 1 eq(3)
225+
B_tilde=np.bmat([[B_,Ac]])
226+
D_tilde=np.bmat([[np.zeros((no_out,no_in)),Cp]])
227+
228+
returnA_,B_tilde,C_,D_tilde,Ap,Bp,Cp
229+
230+
200231
defbalred(sys,orders,method='truncate'):
201232
"""
202233
Balanced reduced order model of sys of a given order.
@@ -222,31 +253,39 @@ def balred(sys, orders, method='truncate'):
222253
Returns
223254
-------
224255
rsys: StateSpace
225-
A reduced order model
256+
A reduced order model or a list of reduced order models if orders is a list
226257
227258
Raises
228259
------
229260
ValueError
230-
* if `method` is not ``'truncate'``
261+
* if `method` is not ``'truncate'`` or ``'matchdc'``
231262
ImportError
232-
if slycot routine ab09ad is not found
263+
if slycot routine ab09ad or ab09bd is not found
264+
265+
ValueError
266+
if there are more unstable modes than any value in orders
233267
234268
Examples
235269
--------
236270
>>> rsys = balred(sys, orders, method='truncate')
237271
238272
"""
239-
ifmethod=='matchdc':
240-
raiseValueError ("MatchDC not yet supported!")
273+
ifmethod!='truncate'andmethod!='matchdc':
274+
raiseValueError("supported methods are 'truncate' or 'matchdc'")
241275
elifmethod=='truncate':
242276
try:
243277
fromslycotimportab09ad
244278
exceptImportError:
245-
raiseControlSlycot("can't find slycot subroutine ab09ad")
246-
else:
247-
raiseValueError("Oops, method is not supported!")
279+
raiseControlSlycot("can't find slycot subroutine ab09ad")
280+
elifmethod=='matchdc':
281+
try:
282+
fromslycotimportab09bd
283+
exceptImportError:
284+
raiseControlSlycot("can't find slycot subroutine ab09bd")
285+
248286

249-
fromscipy.linalgimportschur
287+
fromscipy.linalgimportschur#, cholesky, svd
288+
fromnumpy.linalgimportcholesky,svd
250289
#Check for ss system object, need a utility for this?
251290

252291
#TODO: Check for continous or discrete, only continuous supported right now
@@ -256,6 +295,8 @@ def balred(sys, orders, method='truncate'):
256295
# dico = 'D'
257296
# else:
258297
dico='C'
298+
job='B'# balanced (B) or not (N)
299+
equil='N'# scale (S) or not (N)
259300

260301
rsys= []#empty list for reduced systems
261302

@@ -278,39 +319,26 @@ def balred(sys, orders, method='truncate'):
278319
raiseValueError("System has %i unstable states which is more than ORDER(%i)"% (nn-l,i))
279320

280321
foriinorders:
281-
ifl>0:#handles the stable/unstable decomposition if unstable eigenvalues are found
322+
if(nn-l)>0:#handles the stable/unstable decomposition if unstable eigenvalues are found, nn - l is the number of ustable eigenvalues
282323
#Author: M. Clement (mdclemen@eng.ucsd.edu) 2016
283324
print("Unstable eigenvalues found, performing stable/unstable decomposition")
325+
284326
rorder=i- (nn-l)
285-
As=np.asmatrix(T)
286-
Bs=V.T*sys.B
287-
Cs=sys.C*V
288-
#from ref 1 eq(1) As = [A_ Ac], Bs = [B_], and Cs = [C_ C+]; _ denotes stable subsystem
289-
# [0 A+] [B+]
290-
A_=As[0:l,0:l]
291-
Ac=As[0:l,l::]
292-
Ap=As[l::,l::]
293-
294-
B_=Bs[0:l,:]
295-
Bp=Bs[l::,:]
296-
297-
C_=Cs[:,0:l]
298-
Cp=Cs[:,l::]
299-
#do some more tricky math IAW ref 1 eq(3)
300-
B_tilde=np.bmat([[B_,Ac]])
301-
D_tilde=np.bmat([[np.zeros((rr,mm)),Cp]])
327+
A_,B_tilde,C_,D_tilde,Ap,Bp,Cp=stabsep(T,V,sys,l,mm,rr)
302328

303329
subSys=StateSpace(A_,B_tilde,C_,D_tilde)
304-
305-
job='B'# balanced (B) or not (N)
306-
equil='N'# scale (S) or not (N)
307330
n=np.size(subSys.A,0)
308331
m=np.size(subSys.B,1)
309332
p=np.size(subSys.C,0)
310-
Nr,Ar,Br,Cr,hsv=ab09ad(dico,job,equil,n,m,p,subSys.A,subSys.B,subSys.C,nr=rorder,tol=0.0)
311333

312-
rsubSys=StateSpace(Ar,Br,Cr,np.zeros((p,m)))
334+
ifmethod=='truncate':
335+
Nr,Ar,Br,Cr,hsv=ab09ad(dico,job,equil,n,m,p,subSys.A,subSys.B,subSys.C,nr=rorder,tol=0.0)
336+
rsubSys=StateSpace(Ar,Br,Cr,np.zeros((p,m)))
313337

338+
elifmethod=='matchdc':
339+
Nr,Ar,Br,Cr,Dr,hsv=ab09bd(dico,job,equil,n,m,p,subSys.A,subSys.B,subSys.C,subSys.D,nr=rorder,tol1=0.0,tol2=0.0)
340+
rsubSys=StateSpace(Ar,Br,Cr,Dr)
341+
314342
A_r=rsubSys.A
315343
#IAW ref 1 eq(4) B^{tilde}_r = [B_r, Acr]
316344
B_r=rsubSys.B[:,0:mm]
@@ -324,15 +352,17 @@ def balred(sys, orders, method='truncate'):
324352

325353
rsys.append(StateSpace(Ar,Br,Cr,sys.D))
326354

327-
else:
328-
job='B'# balanced (B) or not (N)
329-
equil='N'# scale (S) or not (N)
355+
else:#stable system branch
330356
n=np.size(sys.A,0)
331357
m=np.size(sys.B,1)
332358
p=np.size(sys.C,0)
333-
Nr,Ar,Br,Cr,hsv=ab09ad(dico,job,equil,n,m,p,sys.A,sys.B,sys.C,nr=i,tol=0.0)
359+
ifmethod=='truncate':
360+
Nr,Ar,Br,Cr,hsv=ab09ad(dico,job,equil,n,m,p,sys.A,sys.B,sys.C,nr=i,tol=0.0)
361+
rsys.append(StateSpace(Ar,Br,Cr,sys.D))
334362

335-
rsys.append(StateSpace(Ar,Br,Cr,sys.D))
363+
elifmethod=='matchdc':
364+
Nr,Ar,Br,Cr,Dr,hsv=ab09bd(dico,job,equil,n,m,p,sys.A,sys.B,sys.C,sys.D,nr=rorder,tol1=0.0,tol2=0.0)
365+
rsys.append(StateSpace(Ar,Br,Cr,Dr))
336366

337367
#if orders was a scalar, just return the single reduced model, not a list
338368
iflen(orders)==1:

‎control/statefbk.py

Lines changed: 38 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -316,7 +316,8 @@ def gram(sys,type):
316316
State-space system to compute Gramian for
317317
type: String
318318
Type of desired computation.
319-
`type` is either 'c' (controllability) or 'o' (observability).
319+
`type` is either 'c' (controllability) or 'o' (observability). To compute the
320+
Cholesky factors of gramians use 'cf' (controllability) or 'of' (observability)
320321
321322
Returns
322323
-------
@@ -327,16 +328,19 @@ def gram(sys,type):
327328
------
328329
ValueError
329330
* if system is not instance of StateSpace class
330-
* if `type` is not 'c'or 'o'
331+
* if `type` is not 'c', 'o', 'cf'or 'of'
331332
* if system is unstable (sys.A has eigenvalues not in left half plane)
332333
333334
ImportError
334-
if slycot routin sb03md cannot be found
335+
if slycot routine sb03md cannot be found
336+
if slycot routine sb03od cannot be found
335337
336338
Examples
337339
--------
338340
>>> Wc = gram(sys,'c')
339341
>>> Wo = gram(sys,'o')
342+
>>> Rc = gram(sys,'cf'), where Wc=Rc'*Rc
343+
>>> Ro = gram(sys,'of'), where Wo=Ro'*Ro
340344
341345
"""
342346

@@ -358,25 +362,42 @@ def gram(sys,type):
358362
foreinD:
359363
ife.real>=0:
360364
raiseValueError("Oops, the system is unstable!")
361-
iftype=='c':
365+
iftype=='c'ortype=='cf':
362366
tra='T'
363-
C=-np.dot(sys.B,sys.B.transpose())
364-
eliftype=='o':
367+
iftype=='c':
368+
C=-np.dot(sys.B,sys.B.transpose())
369+
eliftype=='o'ortype=='of':
365370
tra='N'
366-
C=-np.dot(sys.C.transpose(),sys.C)
371+
iftype=='o':
372+
C=-np.dot(sys.C.transpose(),sys.C)
367373
else:
368374
raiseValueError("Oops, neither observable, nor controllable!")
369375

370376
#Compute Gramian by the Slycot routine sb03md
371377
#make sure Slycot is installed
372-
try:
373-
fromslycotimportsb03md
374-
exceptImportError:
375-
raiseControlSlycot("can't find slycot module 'sb03md'")
376-
n=sys.states
377-
U=np.zeros((n,n))
378-
A=np.array(sys.A)# convert to NumPy array for slycot
379-
X,scale,sep,ferr,w=sb03md(n,C,A,U,dico,job='X',fact='N',trana=tra)
380-
gram=X
381-
returngram
378+
iftype=='c'ortype=='o':
379+
try:
380+
fromslycotimportsb03md
381+
exceptImportError:
382+
raiseControlSlycot("can't find slycot module 'sb03md'")
383+
n=sys.states
384+
U=np.zeros((n,n))
385+
A=np.array(sys.A)# convert to NumPy array for slycot
386+
X,scale,sep,ferr,w=sb03md(n,C,A,U,dico,job='X',fact='N',trana=tra)
387+
gram=X
388+
returngram
389+
eliftype=='cf'ortype=='of':
390+
try:
391+
fromslycotimportsb03od
392+
exceptImportError:
393+
raiseControlSlycot("can't find slycot module 'sb03od'")
394+
n=sys.states
395+
Q=np.zeros((n,n))
396+
A=np.array(sys.A)# convert to NumPy array for slycot
397+
B=np.zeros_like(A)
398+
B[0:sys.B.shape[0],0:sys.B.shape[1]]=sys.B
399+
m=sys.B.shape[0]
400+
X,scale,w=sb03od(n,m,A,Q,B,dico,fact='N',trans=tra)
401+
gram=X
402+
returngram
382403

0 commit comments

Comments
 (0)
[8]ページ先頭
©2009-2025 Movatter.jp