递归最小二乘系统辨识

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['font.sans-serif'] = ['SimHei']
def generate_m_sequence(m):
    register = np.ones(m, dtype=int)
    feedback_poly = [1, 0, 0, 1]  # 例如，x^4 + x + 1
    if len(feedback_poly) != m:
        raise ValueError("Feedback polynomial length must match the register length.")
    m_sequence = []
    for _ in range(2**m - 1):
        m_sequence.append(register[-1])
        new_bit = np.dot(register, feedback_poly) % 2
        register[1:] = register[:-1]
        register[0] = new_bit
    return np.array(m_sequence)

a1_true = 1.5
a2_true = 0.7
b1_true = 1.0
b2_true = 0.5
N = 100

# 生成输入序列 u(k) - 4阶M序列，幅值为5
m = 4
m_sequence = generate_m_sequence(m)
u = 5 * (2 * np.tile(m_sequence, N // len(m_sequence) + 1)[:N] - 1)
V = np.random.normal(0, 1, N)
# 初始化输出序列 z(k)
z = np.zeros(N)
for k in range(2, N):
    z[k] = -a1_true * z[k-1] - a2_true * z[k-2] + b1_true * u[k-1] + b2_true * u[k-2] + V[k]
theta_hat = np.zeros(4)  # 初始参数估计
P = np.eye(4) * 1000  # 初始协方差矩阵，数值较大表示初始不确定性高
lam = 1.0  # 忘记因子
theta_history = []
for k in range(2, N):
    phi_k = np.array([-z[k-1], -z[k-2], u[k-1], u[k-2]])
    K = P @ phi_k / (lam + phi_k.T @ P @ phi_k)
    theta_hat = theta_hat + K * (z[k] - phi_k.T @ theta_hat)
    P = (P - K[:, None] * phi_k @ P) / lam
    theta_history.append(theta_hat.copy())
theta_history = np.array(theta_history)

a1_est, a2_est, b1_est, b2_est = theta_hat

print('最终估计的参数:')
print(f'a1 = {a1_est:.4f}')
print(f'a2 = {a2_est:.4f}')
print(f'b1 = {b1_est:.4f}')
print(f'b2 = {b2_est:.4f}')
plt.figure(figsize=(12, 8))

plt.subplot(4, 1, 1)
plt.plot(theta_history[:, 0], 'b', label='a1估计值')
plt.axhline(y=a1_true, color='r', linestyle='--', label='a1真实值')
plt.legend()
plt.xlabel('迭代次数')
plt.ylabel('a1')

plt.subplot(4, 1, 2)
plt.plot(theta_history[:, 1], 'b', label='a2估计值')
plt.axhline(y=a2_true, color='r', linestyle='--', label='a2真实值')
plt.legend()
plt.xlabel('迭代次数')
plt.ylabel('a2')

plt.subplot(4, 1, 3)
plt.plot(theta_history[:, 2], 'b', label='b1估计值')
plt.axhline(y=b1_true, color='r', linestyle='--', label='b1真实值')
plt.legend()
plt.xlabel('迭代次数')
plt.ylabel('b1')

plt.subplot(4, 1, 4)
plt.plot(theta_history[:, 3], 'b', label='b2估计值')
plt.axhline(y=b2_true, color='r', linestyle='--', label='b2真实值')
plt.legend()
plt.xlabel('迭代次数')
plt.ylabel('b2')

plt.tight_layout()
plt.show()