python實現一般最小二乘系統辨識方法

問題：

對於一個未知參數的系統，往往需要用到系統辨識的方法，例如對於一個單輸入單輸出系統：

　　　　　　　　Z(k)+a₁*Z(k-1)+a₂*Z(k-2)=b₁*U(k-1)+b₂*U(k-2)+V(k)

　　　　　　　　其中：V(k)=c₁*v(k)+c₂*v(k-1)+c₃*v(k-3)

輸入信號采用四階幅值為1的M序列，高斯噪聲v(k)均值為0，方差為0.1。

假設真值為a1=1.6，a2=0.7,b1=1.0,b2=0.4,c1=1.0,c2=1.6,c3=0.7。需要對以上參數進行辨識。

方法：

一般最小二乘法是系統辨識中最簡單的辨識方法之一，其MATLAB實現方法十分簡單，我發現網上一大把，所以我決定“翻譯”一個python版的一般最小二乘辨識方法供讀者參考。

本文用到的庫：

import numpy as np
import matplotlib.pyplot as plt
from operator import xor
from numpy.linalg import inv

 

M序列的產生：

在python中，M序列的產生方法與matlab類似，先產生隨機數，然后采用四階移位寄存器濾波變換，得到我們想要的M序列。

#M序列產生
L=16
#設置M序列周期
#定義初始值
y=np.zeros(L)
u=np.zeros(L)
y1=1
y2=1
y3=1
y4=0
for i in range(0,L):
    x1=xor(y3,y4)  
    x2=y1 
    x3=y2  
    x4=y3 
    y[i]=y4  
    if  y[i]>0.5:
        u[i]=-1  
    else:
        u[i]=1   
    y1=x1
    y2=x2
    y3=x3
    y4=x4     
plt.figure(num=2)
x=np.linspace(0,15,16)
plt.bar(x,u,width=0.1)
plt.title('輸入信號M序列')

高斯分布噪聲：

這里采用的是random庫中的函數，可以看到，我們設置的均值為0，方差為0.1。

 #產生一組16個N(0,1)的高斯分布的隨機噪聲
 mu=0
 sigma=0.1
 samplenum=16
 n=np.random.normal(mu,sigma,samplenum) plt.figure(num=1) plt.plot(n) plt.title("高斯分布隨機噪聲")

 最小二乘辨識程序：

#最小二乘辨識過程
z=np.zeros(16)

for k in range(2,15):   
    z[k]=-1.6*z[k-1]-0.7*z[k-2]+1.0*u[k-1]+0.4*u[k-2]+1.0*n[k]+1.6*n[k-1]+0.7*n[k-2] 
    
plt.figure(num=3) 
plt.bar(x,z,width=0.1)  
plt.title('輸出觀測值')    

H=np.array([[-z[1],-z[0],u[1],u[0]],[-z[2],-z[1],u[2],u[1]],[-z[3],-z[2],u[3],u[2]],[-z[4],-z[3],u[4],u[3]],[-z[5],-z[4],u[5],u[4]],[-z[6],-z[5],u[6],u[5]],[-z[7],-z[6],u[7],u[6]],[-z[8],-z[7],u[8],u[7]],[-z[9],-z[8],u[9],u[8]],[-z[10],-z[9],u[10],u[9]],[-z[11],-z[10],u[11],u[10]],[-z[12],-z[11],u[12],u[11]],[-z[13],-z[12],u[13],u[12]],[-z[14],-z[13],u[14],u[13]]])
Z=np.array([z[2],z[3],z[4],z[5],z[6],z[7],z[8],z[9],z[10],z[11],z[12],z[13],z[14],z[15]]) 

In_1=np.transpose(H)
In_2=np.dot(In_1,H)
In_3=inv(In_2)
In_4=np.dot(In_3,In_1)
c=np.dot(In_4,Z)

分離參數：

#分離參數並顯示
a1=c[0] 
a2=c[1]
b1=c[2]
b2=c[3] 
print("a1的值是：",a1)
print("a2的值是：",a2)
print("b1的值是：",b1)
print("b2的值是：",b2)

注意：

由於在python中plt的庫是不支持中文的，所以要加上這些代碼，保證輸出的圖片標題的中文顯示正常。

#顯示中文字體
plt.rcParams['font.sans-serif']=['SimHei'] #用來正常顯示中文標簽
plt.rcParams['axes.unicode_minus']=False #用來正常顯示負號

結果：

在網上找了半天python畫針狀圖的資料，發現沒有。。所以強行用瘦了的柱狀圖表示針狀圖了。

 

總結：

可以看出Python寫的系統辨識誤差還是有一些的，不過也是受到一般最小二乘參數辨識方法的限制，如果采用遞推最小二乘，增廣最小二乘等方法可能會進一步提高准確性。筆者嘗試過遞推最小二乘，但是與MATLAB相比，其代碼量大大增加，因此在系統辨識方法上，不建議都用Python來寫，MATLAB是個不錯的選擇。當然，喜歡寫python的話，這都不是問題。

代碼全文：

 

# -*- coding: utf-8 -*-
"""
Created on Wed Sep 20 16:11:27 2017
@author: Hangingter
"""
#一般最小二乘辨識

#導入相應科學計算的包
import numpy as np
import matplotlib.pyplot as plt
from operator import xor
from numpy.linalg import inv

#顯示中文字體
plt.rcParams['font.sans-serif']=['SimHei'] #用來正常顯示中文標簽
plt.rcParams['axes.unicode_minus']=False #用來正常顯示負號
#產生一組16個N(0,1)的高斯分布的隨機噪聲
mu=0
sigma=0.1
samplenum=16
n=np.random.normal(mu,sigma,samplenum) 
plt.figure(num=1)
plt.plot(n)
plt.title("高斯分布隨機噪聲")
#M序列產生
L=16
#設置M序列周期
#定義初始值
y=np.zeros(L)
u=np.zeros(L)
y1=1
y2=1
y3=1
y4=0
for i in range(0,L):
    x1=xor(y3,y4)  
    x2=y1 
    x3=y2  
    x4=y3 
    y[i]=y4  
    if  y[i]>0.5:
        u[i]=-1  
    else:
        u[i]=1   
    y1=x1
    y2=x2
    y3=x3
    y4=x4     
plt.figure(num=2)
x=np.linspace(0,15,16)
plt.bar(x,u,width=0.1)
plt.title('輸入信號M序列')
#最小二乘辨識過程
z=np.zeros(16)

for k in range(2,15):   
    z[k]=-1.6*z[k-1]-0.7*z[k-2]+1.0*u[k-1]+0.4*u[k-2]+1.0*n[k]+1.6*n[k-1]+0.7*n[k-2] 
    
plt.figure(num=3) 
plt.bar(x,z,width=0.1)  
plt.title('輸出觀測值')    

H=np.array([[-z[1],-z[0],u[1],u[0]],[-z[2],-z[1],u[2],u[1]],[-z[3],-z[2],u[3],u[2]],[-z[4],-z[3],u[4],u[3]],[-z[5],-z[4],u[5],u[4]],[-z[6],-z[5],u[6],u[5]],[-z[7],-z[6],u[7],u[6]],[-z[8],-z[7],u[8],u[7]],[-z[9],-z[8],u[9],u[8]],[-z[10],-z[9],u[10],u[9]],[-z[11],-z[10],u[11],u[10]],[-z[12],-z[11],u[12],u[11]],[-z[13],-z[12],u[13],u[12]],[-z[14],-z[13],u[14],u[13]]])
Z=np.array([z[2],z[3],z[4],z[5],z[6],z[7],z[8],z[9],z[10],z[11],z[12],z[13],z[14],z[15]]) 

In_1=np.transpose(H)
In_2=np.dot(In_1,H)
In_3=inv(In_2)
In_4=np.dot(In_3,In_1)
c=np.dot(In_4,Z)

#分離參數並顯示
a1=c[0] 
a2=c[1]
b1=c[2]
b2=c[3] 
print("a1的值是：",a1)
print("a2的值是：",a2)
print("b1的值是：",b1)
print("b2的值是：",b2)

 

參考：

MATLAB版的系統辨識一般最小二乘方法：

http://blog.csdn.net/sinat_20265495/article/details/51426537