递推最小二乘参数估计RLS

系统辨识与自适应控制MATLAB仿真 修订版  

仿真实例 2.6 递推最小二乘法估计

import numpy as np 
import matplotlib.pyplot as plt 
from mxulie import M_sequences

if __name__ == '__main__':
	L = 20 #序列长度
	Y = np.zeros(L)
	phi = np.zeros((L,4))
	[M,IM]=M_sequences(L)
	xi = np.sqrt(0) * np.random.randn(L,1)
	
	y1 = y2 =0
	u4 = u3 = u2 =u1 =0
	theta = np.array([1.5,-0.7,1,0.5])
	P = 1e6 * np.eye(4)
	theta1 = np.zeros((L,4))
	theta1_1 = np.zeros(4)

	for i in np.arange(L):
		phi[i,:] = np.array([y1 , y2 ,u3 ,u4])		
		#Y[i] = 1.5*y1 -0.7*y2 + u3 + 0.5*u4 + xi[i]
		Y[i] = np.dot(theta,phi[i,:]) + xi[i]


		# RLS       
		K = np.dot(P,phi[i,:])/(1+np.dot(np.dot(phi[i,:],P),phi[i,:]))
		theta1[i,:] = theta1_1 + K*(Y[i]-np.dot(phi[i,:],theta1_1))
		P = np.dot(np.eye(4)-phi[i,:]*K.reshape((-1,1)),P)


		# 数据更新
		theta1_1 = theta1[i,:]
		y2 = y1
		y1 = Y[i]

		u4 = u3
		u3 = u2
		u2 = u1
		u1 = IM[i]

	
	plt.subplot(3,1,1)
	plt.title('输入-逆M序列')
	plt.step(np.arange(L),IM)

	plt.subplot(3,1,2)
	plt.title('输出-Y')
	plt.plot(np.arange(L),Y)

	plt.subplot(3,1,3)
	plt.title('系数Theta')   
	plt.plot(theta1)

	plt.subplots_adjust(hspace = 0.5)
	plt.show()