Python实现系统辨识LS；一般最小二乘；递推最小二乘；

`#M序列产生**

**```#M序列产生**
import numpy as np
import matplotlib.pyplot as plt
from operator import xor
from numpy.linalg import inv
#L=int(input('序列周期'))
L=500
y=np.zeros(L)
u=np.zeros(L)
#500默认
Q=9
Y=np.zeros(L)
X=np.zeros(L)

for i in range(Q-1):
    Y[i]=1

for i in range(0,L):
    X[0]=xor(int(Y[Q-2]),int(Y[Q-1])) 
    for j in range(1,Q):
        X[j]=Y[j-1]
    y[i]=Y[Q-1]  
    if  y[i]>0.5:
        u[i]=-1  
    else:
        u[i]=1   
    for j in range(Q):
        Y[j]=X[j] 
plt.figure(num=2)
x=np.linspace(0,15,L)
print(u)
plt.bar(x,u,width=0.1)
plt.title('12')
a = np.arange(0,L,1)
for i in range(L):
    a[i]=u[i]

np.savetxt('C:\\Users\\****\\Desktop\\DU.txt',a,fmt='%d')

M部分结果

#产生一组500个N(0,1)的高斯分布的随机噪声

import numpy as np
n=np.random.normal(0,0.1,500) 
np.savetxt('C:\\Users\\****\\N.txt',n,fmt='%f')

#观测数据

import numpy as np
z=np.zeros(500)
u=np.loadtxt('C:\\Users\\****\\DU.txt',dtype=np.double)
n=np.loadtxt('C:\\Users\\****\\N.txt',dtype=np.float)

for k in range(2,499):   
    z[k]=1.5*z[k-1]-0.7*z[k-2]+1.0*u[k-1]+0.5*u[k-2]+1.0*n[k]+1.6*n[k-1]+0.7*n[k-2] 
np.savetxt('C:\\Users\\****\\Z.txt',z,fmt='%f')

#一般LS 辨识参数

import numpy as np
import matplotlib.pyplot as plt
from operator import xor
from numpy.linalg import inv
import numpy as np
Z=np.loadtxt('C:\\Users\\涛片天\\Desktop\\Z.txt',dtype=np.float32)
U=np.loadtxt('C:\\Users\\涛片天\\Desktop\\DU.txt',dtype=np.double)

H=np.zeros((498,4))
for i in range(0,498):
    H[i,0]=Z[i]
    H[i,1]=Z[i-1]
    H[i,2]=U[i]
    H[i,3]=U[i-1]

Q=Z[2:]
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,Q)
a1=c[0] 
a2=c[1]
b1=c[2]
b2=c[3] 
print("a1的值是：",a1)
print("a2的值是：",a2)
print("b1的值是：",b1)
print("b2的值是：",b2)

``

#递推LS

import numpy as np
import matplotlib.pyplot as plt
from operator import xor
from numpy.linalg import inv
import numpy as np
Z=np.loadtxt('C:\\Users\\***\\Z.txt',dtype=np.float32)
U=np.loadtxt('C:\\Users\\***\\DU.txt',dtype=np.double)
c0=np.zeros((4,1)) 
h1=np.zeros((4,1))   
c0=c0+0.1
p0=np.eye(4)
E=np.eye(4)
p0=1000*p0
c=np.zeros((4,499))
lamt=1
for k in range(2,499):
    h1[0,0]=Z[k-1]
    h1[1,0]=Z[k-2]
    h1[2,0]=U[k-1]
    h1[3,0]=U[k-2]
    m=inv((h1.T).dot(p0).dot(h1)+lamt)
    k1=p0.dot(h1)
    k1=k1.dot(m)
    new=Z[k]-h1.T.dot(c0)
    c1=c0+k1.dot(new)
    p1=(E-k1.dot(h1.T)).dot(p0)
    c[0:,k:k+1]=c1
    c0=c1
    p0=p1
res=c1
print(res)
np.savetxt('C:\\Users\\***\\res.txt',res,fmt='%f')