2012-10-23 177 views
0

我有一個卡爾曼濾波器的實現,我正在做矩陣運算。在某些時候,我應該減去兩個1x1矩陣。我有一個錯誤,我不知道它是從哪裏來的。 在python中做矩陣運算的最好方法是什麼?python/numpy中的矩陣相減

import numpy as np 
import pylab as pl 
import scipy as Sci 
import scipy.linalg as linalg 

class GetPos(object): 
    def __init__(self): 
     self.Posp = 0 
     self.Velp = 80 
     self.z = np.matrix(0) 
    def __repr__(self): 
     return "from GetPos.__repr__ z=%s" % (self.z) 
    def __call__(self): 
     self.dt = 0.1 
     self.w = 0 + 10*np.random.random() 
     self.v = 0 + 10*np.random.random()    
     self.z = self.Posp + self.Velp*self.dt + self.v 
     self.Posp = self.z - self.v 
     self.Velp = 80 + self.w 
     print 'from GetPos.__call__ z = %s' % self.z 
     return self.z 

class DvKalman(object): 
    def __init__(self): 
     self.dt = .1 
     self.A = np.matrix([[1., self.dt],[0,1]]) 
     self.H = np.matrix([1., 0]) 
     self.Q = np.matrix([[1,0.],[0,3]]) 
     self.R = np.matrix(10) 
     self.x = np.matrix([0,20]).T 
     self.P = np.matrix(5*np.eye(2)) 
     #print 'P matrix \n%s' % self.P 
     self.firstRun = 0 
    def __call__(self, z): 
     self.z = z 
     print 'from DvKalman.__call__ slef.z = %s and z = %s' % (self.z,z) 
     self.xp = self.A * self.x 
     self.Pp = self.A*self.P*self.A.T + self.Q 
     self.K = self.Pp * self.H.T * linalg.inv(np.absolute(self.H*self.Pp*self.H.T + self.R)); 
     print 'from DvKalman.__call__ z=%s, \npreviouse x=\n%s \nH = \n%s \nand P=\n%s \nand xp=\n%s,\n Pp = \n%s,\n K=\n%s' % (self.z,self.x,self.H, self.P,self.xp,self.Pp,self.K) 
     newM1 = self.H*self.xp  
     print 'This is self.H*self.xp %s and this is self.z = %s' % (newM1, self.z) 
     newM2 = self.z - self.H*self.xp 
     print 'This should give simple substruction %s' % newM2     
     self.x = self.xp + self.K*(self.z - self.H*self.xp) 
     self.P = self.Pp - self.K*self.H*self.Pp 
     print 'new values x=%s and P=%s' % (self.x,self.P) 
     return (self.x) 
def TestDvKalman(): 
    Nsamples = np.arange(0,10,.1) 

    kal = DvKalman() 
    #print type(kal) 
    Xsaved = [] 
    Zsaved = [] 

    for i in range(len(Nsamples)): 
     z = GetPos() 
     print z 
     print 'from TestDvKalman zpos = %s' % z 
     Zsaved.append(z) 
     [position, velocity] = kal(z) 
     print position, velocity 
     Xsaved.append([position, velocity]) 
    print Zsaved 
    print Xsaved 
# f1 = pl.subplot(121) 
# f1 = pl.plot(Xsaved, 'x-',label = 'Xsaved') 
# f1 = pl.legend() 
#  
# f2 = pl.subplot(122) 
# f2 = pl.title('Kalman Velocity') 
# f2 = pl.plot(Zsaved, 'o-', color = 'brown',label = 'Zsaved') 
# f2 = pl.legend() 
#  
# pl.show() 

if __name__ == '__main__': 
    TestDvKalman() 

我已經添加了幾行print跟蹤和調試代碼,我增加了新的變數newM這不會是在代碼中。矩陣打印正確This is self.H*self.xp [[ 2.]] and this is self.z = from GetPos.__repr__ z=[[0]]這兩個矩陣都是1x1,但我仍然收到一個錯誤,不知道爲什麼。錯誤是:

newM2 = self.z - self.H*self.xp 
TypeError: unsupported operand type(s) for -: 'GetPos' and 'matrix' 

我懷疑我搞亂了類型的地方,但不知道在哪裏以及如何糾正它。你能指出我的錯誤在哪裏,以及如何構建這樣的代碼以避免將來出現類似的錯誤?

回答

1

通過newM2 = self.z() - self.H*self.xp更換newM2 = self.z - self.H*self.xp

該程序應該與那個(但我不能確認它是否會想要你想要)

2

您正在將GetPos實例傳遞給DvKalman __call__方法。所以你試圖減去一個GetPos實例和一個矩陣。不是矩陣和矩陣。

+0

謝謝。這是正確的,但我應該怎麼做比。每次我給'DvKalman'打電話時,我都想從'GetPos'得到新的變量'z' – tomasz74

2

TestDvKalman,這條線

z = GetPos() 

z到的GetPos一個實例。您在此行中使用這個作爲參數傳遞給kal

[position, velocity] = kal(z) 

所以DvKalman__call__方法給出的GetPos一個實例,其中保存爲self.z。這導致錯誤在這行:

newM2 = self.z - self.H*self.xp 
+0

謝謝。我打算做這樣的事情,但我想''=''[m]]'從'GetPos .__ call__'返回'm'我明白這是它的工作原理。我錯過了什麼? – tomasz74

+1

GetPos()創建一個GetPos的實例。然後您必須*調用GetPos()返回的對象。無論你在創建z時還是在DvKalman的__call__方法中都這樣做,取決於你希望如何工作。 –