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

Commit8b337f0

Browse files
committed
Implement okid
1 parent9258a0c commit8b337f0

File tree

2 files changed

+38
-44
lines changed

2 files changed

+38
-44
lines changed

‎control/modelsimp.py‎

Lines changed: 33 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -561,8 +561,8 @@ def markov(Y, U, m=None, transpose=False):
561561
defokid(*args,m=None,transpose=False,dt=True,truncate=False):
562562
"""okid(Y, U, [, m])
563563
564-
Calculate the first `m` Markov parameters [D CB CAB ...]
565-
from data
564+
Calculate the first `m+1` Markov parameters [D CB CAB ...]
565+
from data.
566566
567567
This function computes the Markov parameters for a discrete time system
568568
@@ -578,12 +578,12 @@ def okid(*args, m=None, transpose=False, dt=True, truncate=False):
578578
579579
The function can be called with either 1, 2 or 3 arguments:
580580
581-
* ``H = okid(response)``
582-
* ``H = okid(respnose, m)``
581+
* ``H = okid(data)``
582+
* ``H = okid(data, m)``
583583
* ``H = okid(Y, U)``
584584
* ``H = okid(Y, U, m)``
585585
586-
where `response` is an `TimeResponseData` object, and `Y`, `U`, are 1D or 2D
586+
where `data` is an `TimeResponseData` object, and `Y`, `U`, are 1D or 2D
587587
array and m is an integer.
588588
589589
Parameters
@@ -601,10 +601,10 @@ def okid(*args, m=None, transpose=False, dt=True, truncate=False):
601601
m : int, optional
602602
Number of Markov parameters to output. Defaults to len(U).
603603
dt : True of float, optional
604-
True indicates discrete time with unspecified sampling time,
605-
positivenumber is discrete time with specified sampling time.It
606-
can be used to scale themarkov parameters in order to match the
607-
impulse response ofthis library. Default is True.
604+
True indicates discrete time with unspecified sampling time and a
605+
positivefloat is discrete time withthespecified sampling time.
606+
Itcan be used to scale theMarkov parameters in order to match
607+
the unit-areaimpulse response ofpython-control. Default is True
608608
truncate : bool, optional
609609
Do not use first m equation for least least squares. Default is False.
610610
transpose : bool, optional
@@ -615,7 +615,7 @@ def okid(*args, m=None, transpose=False, dt=True, truncate=False):
615615
Returns
616616
-------
617617
H : ndarray
618-
First m Markov parameters, [D CB CAB ...]
618+
First m Markov parameters, [D CB CAB ...].
619619
620620
References
621621
----------
@@ -683,29 +683,17 @@ def okid(*args, m=None, transpose=False, dt=True, truncate=False):
683683

684684
# the algorithm - Construct a matrix of control inputs to invert
685685
#
686-
# (q,l) = (q,p*m) @ (p*m,l)
687-
# YY.T =H @UU.T
686+
# (q,l) = (q,(p+q)*m+p) @ ((p+q)*m+p,l)
687+
# YY.T =Ybar @VV.T
688688
#
689689
# This algorithm sets up the following problem and solves it for
690690
# the Markov parameters
691691
#
692-
# (l,q) = (l,p*m) @ (p*m,q)
693-
# YY = UU @ H.T
694-
#
695-
# [ y(0) ] [ u(0) 0 0 ] [ D ]
696-
# [ y(1) ] [ u(1) u(0) 0 ] [ C B ]
697-
# [ y(2) ] = [ u(2) u(1) u(0) ] [ C A B ]
698-
# [ : ] [ : : : : ] [ : ]
699-
# [ y(l-1) ] [ u(l-1) u(l-2) u(l-3) ... u(l-m) ] [ C A^{m-2} B ]
692+
# (l,q) = (l,(p+q)*m+p) @ ((p+q)*m+p,q)
693+
# YY = VV @ Ybar.T
700694
#
701695
# truncated version t=m, do not use first m equation
702696
#
703-
# [ y(t) ] [ u(t) u(t-1) u(t-2) u(t-m) ] [ D ]
704-
# [ y(t+1) ] [ u(t+1) u(t) u(t-1) u(t-m+1)] [ C B ]
705-
# [ y(t+2) ] = [ u(t+2) u(t+1) u(t) u(t-m+2)] [ C B ]
706-
# [ : ] [ : : : : ] [ : ]
707-
# [ y(l-1) ] [ u(l-1) u(l-2) u(l-3) ... u(l-m) ] [ C A^{m-2} B ]
708-
#
709697
# Note: This algorithm assumes C A^{j} B = 0
710698
# for j > m-2. See equation (3) in
711699
#
@@ -716,42 +704,46 @@ def okid(*args, m=None, transpose=False, dt=True, truncate=False):
716704
#
717705

718706
Vmat=np.concatenate((Umat,Ymat),axis=0)
719-
print(Vmat.shape)
720707

721-
VVT=np.zeros(((q+p)*m+p,l))
708+
VVT=np.zeros(((p+q)*m+p,l))
722709

723710
VVT[:p,:]=Umat
724711
foriinrange(m):
725712
# Shift previous column down and keep zeros at the top
726-
VVT[i*(q+p)+p:(i+1)*(q+p)+p,i+1:]=Vmat[:,:l-i-1]
713+
VVT[(p+q)*i+p:(p+q)*(i+1)+p,i+1:]=Vmat[:,:l-i-1]
727714

728715
YY=Ymat[:,t:].T
729716
VV=VVT[:,t:].T
730717

731-
# Solve for the Markov parameters from YY = VV @ Ybar.T
718+
# Solve for theobserverMarkov parameters from YY = VV @ Ybar.T
732719
YbarT,_,_,_=np.linalg.lstsq(VV,YY,rcond=None)
733-
Ybar=YbarT.T/dt# scaling
720+
Ybar=YbarT.T
734721

722+
# Equation 11, Page 9
735723
D=Ybar[:,:p]
736-
print(D.shape)
737724

738-
Ybar_r=Ybar[:,p:].reshape(q,m,q+p)# output, time*input -> output, time,input
739-
Ybar_r=Ybar_r.transpose(0,2,1)# output,input, time
725+
Ybar_r=Ybar[:,p:].reshape(q,m,p+q)# output, time*v_input -> output, time,v_input
726+
Ybar_r=Ybar_r.transpose(0,2,1)# output,v_input, time
740727

741-
Ybar1=Ybar_r[:,:p,:]#from input
742-
Ybar2=Ybar_r[:,p:,:]#from output
728+
Ybar1=Ybar_r[:,:p,:]#select observer Markov parameters generated by input
729+
Ybar2=Ybar_r[:,p:,:]#select observer Markov parameters generated by output
743730

731+
# Using recursive substitution to solve for Markov parameters
744732
Y=np.zeros((q,p,m))
745-
Y[:,:,0]=Ybar1[:,:,0]+Ybar1[:,:,0]@D
733+
# Equation 12, Page 9
734+
Y[:,:,0]=Ybar1[:,:,0]+Ybar2[:,:,0]@D
735+
736+
# Equation 15, Page 10
746737
forkinrange(1,m):
747-
Y[:,:,i]=Ybar1[:,:,i]+Ybar1[:,:,i]@D
738+
Y[:,:,k]=Ybar1[:,:,k]+Ybar2[:,:,k]@D
748739
foriinrange(k-1):
749740
Y[:,:,k]+=Ybar2[:,:,i]@Y[:,:,k-i-1]
750741

751-
752-
H=np.zeros((q,p,m))
742+
# Equation 11, Page 9
743+
H=np.zeros((q,p,m+1))
753744
H[:,:,0]=D
754-
H[:,:,1:]=Y[:,:,:-1]
745+
H[:,:,1:]=Y[:,:,:]
746+
H=H/dt# scaling
755747

756748
# Return the first m Markov parameters
757749
returnHifnottransposeelsenp.transpose(H)

‎examples/okid_msd.py‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -67,7 +67,7 @@ def create_impulse_response(H, time, transpose, dt):
6767

6868

6969
xixo_list= ["SISO","SIMO","MISO","MIMO"]
70-
xixo=xixo_list[0]# choose a system for estimation
70+
xixo=xixo_list[3]# choose a system for estimation
7171
matchxixo:
7272
case"SISO":
7373
sys=ct.StateSpace(A,B[:,0],C[0,:],D[0,0])
@@ -90,8 +90,10 @@ def create_impulse_response(H, time, transpose, dt):
9090
response.plot()
9191
plt.show()
9292

93-
m=200
94-
ir_true=ct.impulse_response(sysd,T=dt*m)
93+
m=100
94+
ir_true=ct.impulse_response(sysd,T=t[:m+1])
95+
H_true=ir_true.outputs
96+
print(H_true.shape)
9597

9698
H_est=ct.okid(response,m,dt=dt)
9799
print(H_est.shape)

0 commit comments

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