自动控制 一阶线性系统参数测定

本文已参与「新人创作礼」活动，一起开启掘金创作之路。

想做一个PID温控电烙铁，首先需要测定系统参数，便于以后的仿真和选择PID参数 设电烙铁发热器的发热功率为$x(t)$（输入），比外界温度高为$y(t)$（输出），热导为$G$，热容量为$C$，可列出以下微分方程： $C\frac{dy(t)}{dt}+Gy(t)=x(t)$ 时间常数$\tau=\frac{C}{G}$ 该微分方程的通解$y(t)=y_0+\frac{1}{C}e^{-t/\tau}\int_{-\infty}^{t}x(t)e^{t/\tau}dt$ 先来看看仿真结果

# 参数
C = 20
G = 1
t = np.arange(0, 600, 0.1)
Pf = lambda t: 1 + np.sin(t**2/5000)
P = Pf(t) #发热功率必须为正数

# 生成数据，加入噪声
#pn = P.mean()/50
#P += np.random.normal(0, pn, P.shape)  # 实际功率噪声,odeint无效
# 使用通解进行计算，但数值太大，可能会浮点溢出，不建议使用
T = 1/C*np.exp(-G/C*t)*np.cumsum(P*np.exp(G/C*t))*(t[1]-t[0])
# 使用scipy求解ODE，不能仿真实际功率噪声
#T = odeint(lambda T,t_:(Pf(t_)-T*G)/C, 0, t)[:, 0]
pn = P.mean()/50
P += np.random.normal(0, pn, P.shape)  # 功率测量噪声
tn = T.mean()/50
T += np.random.normal(0, tn, T.shape)  # 温度测量噪声

# 温度测量比功率测量落后(>0)或超前(<0)多少帧
dt = 0
if dt > 0:
    P, T, t = P[:-dt], T[dt:], t[:-dt]
if dt < 0:
    P, T, t = P[-dt:], T[:dt], t[:dt]

# 绘图
plt.plot(t, P, c='r')
plt.twinx()
plt.plot(t, T, c='b')
plt.show()

傅里叶变换一下，看看频谱

Pw = np.fft.rfft(P)
Tw = np.fft.rfft(T)
PwTw = Pw/Tw
w = np.linspace(0, 2*np.pi/(2*(t[1]-t[0])), Pw.shape[0])

n = 30
plt.plot(w[:n], np.log(np.abs(Pw[:n])), c='r')
plt.twinx()
plt.plot(w[:n], np.log(np.abs(Tw[:n])), c='b')
plt.show()

经过傅里叶变换，变为 $j\omega CY(\omega)+GY(\omega)=X(\omega)$ $(j\omega C+G)Y(\omega)=X(\omega)$ 实部$G *Re[Y(\omega)]=Re[X(\omega)]$ 虚部$\omega C*Im[Y(\omega)]=Im[X(\omega)]$ 实际测量到的数据是有噪声的 $H^{-1}(\omega)=\frac{X(\omega)}{Y(\omega)}=j\omega C+G$ $G=Re[H^{-1}(\omega)]$ $C=\frac{1}{\omega}Im[H^{-1}(\omega)]$

n = 30
plt.plot(w[:n], PwTw.real[:n], c='r')
plt.twinx()
plt.plot(w[:n], PwTw.imag[:n], c='b')
plt.twiny()
plt.xlim(-1, n+1)  # 数组索引
plt.show()

只要X(s)输入啁啾信号，这样可以在实验中测定$C$和$G$ 另外，实际测量到的数据是有噪声的，为了提高分辨率，需要进行最小二乘法 $Q(G, U)=\sum\left|X_i-(j\omega_iC+G)Y_i\right|^2=\sum\left|\varepsilon_i\right|^2=\sum[Re^2(\varepsilon_i)+Im^2(\varepsilon_i)]$