小西作业1_third order plant(SPM)
Consider the third order plant
y
=
G
(
x
)
u
y=G(x)u
y=G(x)u
where
G
(
s
)
=
b
2
s
2
+
b
1
s
+
b
0
s
3
+
a
2
s
2
+
a
1
s
+
a
0
G(s)=\dfrac{b_2s^2+b_1s+b_0}{s^3+a_2s^2+a_1s+a_0}
G(s)=s3+a2s2+a1s+a0b2s2+b1s+b0
(a) For the case where all parameters are unknown, obtain a static linear parametric model (SPM) of the plant. Assume that
u
u
u and
y
y
y are available for measurement, but their derivatives are not.
(b) For the case where it is known that
a
0
=
b
2
=
2
a_0=b_2=2
a0=b2=2, obtain a parametric model for the plant in terms of
θ
∗
=
[
b
1
,
b
0
,
a
2
,
a
1
]
T
\theta^*=[b_1,b_0,a_2,a_1]^T
θ∗=[b1,b0,a2,a1]T.
(c) For the case where it is known that
b
0
=
b
1
=
0
,
b
2
=
3
b_0=b_1=0, b_2=3
b0=b1=0,b2=3 obtain a parametric model in terms of
θ
∗
=
[
a
2
,
a
1
,
a
0
]
T
\theta^*=[a_2,a_1,a_0]^T
θ∗=[a2,a1,a0]T. Design and simulate in Matlab/Simulink a gradient based algorithm to generate on line esrimates of
θ
∗
\theta^*
θ∗. For the actual system you can take
a
0
=
b
2
=
2
,
b
0
=
b
1
=
0
,
b
2
=
3
,
a
1
=
a
2
=
2
a_0=b_2=2,b_0=b_1=0,b_2=3,a_1=a_2=2
a0=b2=2,b0=b1=0,b2=3,a1=a2=2. First use
u
(
t
)
=
1
u(t)=1
u(t)=1, then use
u
(
t
)
=
sin
(
0.2
t
)
+
sin
(
t
)
u(t)=\sin(0.2t)+\sin(t)
u(t)=sin(0.2t)+sin(t). Compare the different results.
Solution:
(a) both sides time
s
3
+
a
2
s
2
+
a
1
s
+
a
0
s^3 + a_2s^2+a_1s+a_0
s3+a2s2+a1s+a0, then
(
s
3
+
a
2
s
2
+
a
1
s
+
a
0
)
y
=
(
b
2
s
2
+
b
1
s
+
b
0
)
u
(s^3+a_2s^2+a_1s+a_0)y=(b_2s^2+b_1s+b_0)u
(s3+a2s2+a1s+a0)y=(b2s2+b1s+b0)uso we get the solution
y
′
′
′
+
a
2
y
′
′
+
a
1
y
′
+
a
0
y
=
b
2
u
′
′
+
b
1
u
′
+
b
0
u
y'''+a_2y''+a_1y'+a_0y=b_2u''+b_1u'+b_0u
y′′′+a2y′′+a1y′+a0y=b2u′′+b1u′+b0u
then we can obtain a static linear parametric model (SPM) of the plant
y
′
′
′
=
θ
∗
T
ϕ
y'''={\theta^*}^T\phi
y′′′=θ∗Tϕ
where
θ
∗
=
[
a
2
a
1
a
0
b
2
b
1
b
0
]
,
ϕ
=
[
−
y
′
′
−
y
′
−
y
u
′
′
u
′
u
]
{\theta^*}=\begin{bmatrix}a_2\\ a_1\\ a_0\\ b_2\\ b_1\\ b_0 \end{bmatrix},\phi=\begin{bmatrix}-y''\\ -y'\\ -y\\ u''\\ u'\\ u \end{bmatrix}
θ∗=
a2a1a0b2b1b0
,ϕ=
−y′′−y′−yu′′u′u
since
u
u
u and
y
y
y are available for measurement, but their derivatives are not, then we get
a
0
y
=
b
0
u
a_0y=b_0u
a0y=b0u
So the static linear parametric mode of the plant is
y
=
b
0
a
0
u
y=\dfrac{b_0}{a_0}u
y=a0b0u
(b) when
a
0
=
b
0
=
2
a_0=b_0=2
a0=b0=2
G
(
s
)
=
2
s
2
+
b
1
s
+
b
0
s
3
+
a
2
s
2
+
a
1
s
+
2
G(s)=\dfrac{2s^2+b_1s+b_0}{s^3+a_2s^2+a_1s+2}
G(s)=s3+a2s2+a1s+22s2+b1s+b0then
y
=
2
s
2
+
b
1
s
+
b
0
s
3
+
a
2
s
2
+
a
1
s
+
2
u
y=\dfrac{2s^2+b_1s+b_0}{s^3+a_2s^2+a_1s+2}u
y=s3+a2s2+a1s+22s2+b1s+b0uwe get
(
s
3
+
a
2
s
2
+
a
1
s
+
2
)
y
=
(
2
s
2
+
b
1
s
+
b
0
)
u
(s^3+a_2s^2+a_1s+2)y=(2s^2+b_1s+b_0)u
(s3+a2s2+a1s+2)y=(2s2+b1s+b0)uthen
s
3
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
y
⏟
z
+
a
2
s
2
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
y
+
a
1
s
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
y
+
2
1
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
y
−
2
s
2
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
u
−
b
1
s
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
u
−
b
0
1
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
u
=
0
\begin{equation*}\begin{split}&\underbrace{\dfrac{s^3}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}y}_{z}+a_2\dfrac{s^2}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}y+a_1\dfrac{s}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}y+2\dfrac{1}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}y\\-&2\dfrac{s^2}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}u-b_1\dfrac{s}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}u-b_0\dfrac{1}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}u=0\end{split}\end{equation*}
−z
s3+λ2s2+λ1s+λ0s3y+a2s3+λ2s2+λ1s+λ0s2y+a1s3+λ2s2+λ1s+λ0sy+2s3+λ2s2+λ1s+λ01y2s3+λ2s2+λ1s+λ0s2u−b1s3+λ2s2+λ1s+λ0su−b0s3+λ2s2+λ1s+λ01u=0
then we can obtain a static linear parametric model (SPM) of the plant
z
=
θ
∗
T
ϕ
+
2
1
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
y
−
2
s
2
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
u
z={\theta^*}^T\phi+2\dfrac{1}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}y-2\dfrac{s^2}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}u
z=θ∗Tϕ+2s3+λ2s2+λ1s+λ01y−2s3+λ2s2+λ1s+λ0s2u
where
θ
∗
=
[
a
2
a
1
b
1
b
0
]
,
ϕ
=
[
−
s
2
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
y
−
s
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
y
s
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
u
1
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
u
]
\theta^*=\begin{bmatrix}a_2\\ a_1\\ b_1\\ b_0 \end{bmatrix},\phi=\begin{bmatrix}-\frac{s^2}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}y\\ -\frac{s}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}y\\ \frac{s}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}u\\ \frac{1}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}u \end{bmatrix}
θ∗=
a2a1b1b0
,ϕ=
−s3+λ2s2+λ1s+λ0s2y−s3+λ2s2+λ1s+λ0sys3+λ2s2+λ1s+λ0sus3+λ2s2+λ1s+λ01u
since
u
u
u and
y
y
y are available for measurement, but their derivatives are not, then we get
2
y
=
b
0
u
2y=b_0u
2y=b0u
(c) when
b
0
=
b
1
=
0
,
b
2
=
3
b_0=b_1=0,b_2=3
b0=b1=0,b2=3
G
(
s
)
=
3
s
2
s
3
+
a
2
s
2
+
a
1
s
+
a
0
G(s)=\dfrac{3s^2}{s^3+a_2s^2+a_1s+a_0}
G(s)=s3+a2s2+a1s+a03s2
then
y
=
3
s
2
s
3
+
a
2
s
2
+
a
1
s
+
a
0
u
y=\dfrac{3s^2}{s^3+a_2s^2+a_1s+a_0}u
y=s3+a2s2+a1s+a03s2uso we get
(
s
3
+
a
2
s
2
+
a
1
s
+
a
0
)
y
=
3
s
2
u
(s^3+a_2s^2+a_1s+a_0)y=3s^2u
(s3+a2s2+a1s+a0)y=3s2uthen
s
3
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
y
⏟
z
+
a
2
s
2
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
y
+
a
1
s
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
y
+
a
0
1
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
y
−
3
s
2
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
u
=
0
\begin{equation*} \begin{split}& \underbrace{\dfrac{s^3}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}y}_z+a_2\dfrac{s^2}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}y\\+&a_1\dfrac{s}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}y+a_0\dfrac{1}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}y-3\dfrac{s^2}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}u=0 \end{split} \end{equation*}
+z
s3+λ2s2+λ1s+λ0s3y+a2s3+λ2s2+λ1s+λ0s2ya1s3+λ2s2+λ1s+λ0sy+a0s3+λ2s2+λ1s+λ01y−3s3+λ2s2+λ1s+λ0s2u=0
then we can obtain a parametric model
z
=
θ
∗
T
ϕ
−
3
s
2
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
u
z={\theta^*}^T\phi-3\dfrac{s^2}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}u
z=θ∗Tϕ−3s3+λ2s2+λ1s+λ0s2u
where
θ
∗
=
[
a
2
a
1
a
0
]
,
ϕ
=
[
−
s
2
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
y
−
s
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
y
−
1
s
3
+
λ
2
s
2
+
λ
1
s
+
λ
0
y
]
\theta^*=\begin{bmatrix}a_2\\ a_1\\ a_0 \end{bmatrix},\phi=\begin{bmatrix}-\frac{s^2}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}y\\ -\frac{s}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}y\\ -\frac{1}{s^3+\lambda_2s^2+\lambda_1s+\lambda_0}y \end{bmatrix}
θ∗=
a2a1a0
,ϕ=
−s3+λ2s2+λ1s+λ0s2y−s3+λ2s2+λ1s+λ0sy−s3+λ2s2+λ1s+λ01y
Design and simulate in Matlab a gradient based algorithm to generate on line estimates of θ ∗ \theta^* θ∗. For the actual system we take a 0 = 2 , b 0 = b 1 = 0 , b 2 = 3 , a 1 = a 2 = 2 a_0 =2, b_0 = b_1 = 0, b_2 = 3,a_1 = a_2 = 2 a0=2,b0=b1=0,b2=3,a1=a2=2. First use u ( t ) = 1 u(t) = 1 u(t)=1, then use u ( t ) = sin ( 0.2 t ) + sin ( t ) u(t) = \sin (0.2t) + \sin (t) u(t)=sin(0.2t)+sin(t). Compare the different results.
// An highlighted block
clear,clc
% 定义系统真实参数
a0_true = 2; a1_true = 2; a2_true = 2;
% 初始化参数估计值
theta_hat = [1; 1; 1];
alpha = 0.01;% 定义学习率
epsilon = 1e-6; % 小扰动值
% 定义模拟时间步长和总时间
dt = 0.1; T = 100; t = 0:dt:T;
% 定义输入信号
% u = ones(size(t));%u(t)=1
u = sin(0.2*t) + sin(t);%
u = u(:); % 将输入信号转换为列向量
% 初始化输出估计值和误差
y_hat = zeros(size(t));
e = zeros(size(t));
% 根据实际参数计算整个时间序列的估计输出
G_true = tf([3 0 0], [1 2 2 2]);
y_true = lsim(G_true, u, t); % 一次性计算整个时间序列的估计输出
for k = 1:1
% 根据当前估计参数计算整个时间序列的估计输出
G_hat = tf([3 0 0], [1 theta_hat(1) theta_hat(2) theta_hat(3)]);
y_hat = lsim(G_hat, u, t); % 一次性计算整个时间序列的估计输出
% 计算整个时间序列的误差
e = y_true - y_hat; % 这里假设实际输出为y_true,根据系统传递函数关系
% 计算梯度(使用有限差分法)
J_theta2_original = (e.^2)/2; % 损失函数,误差的平方
grad_J = zeros(size(theta_hat));
for i = 1:length(theta_hat)
theta_hat_perturbed = theta_hat;
theta_hat_perturbed(i) = theta_hat_perturbed(i) + epsilon;
G_hat_perturbed = tf([3 0 0], [1 theta_hat_perturbed(1) theta_hat_perturbed(2) theta_hat_perturbed(3)]);
y_hat_perturbed = lsim(G_hat_perturbed, u, t);
e_perturbed = y_true - y_hat_perturbed;
J_theta2_perturbed = (e_perturbed.^2)/2;
grad_J(i) = mean((J_theta2_perturbed - J_theta2_original) / epsilon); % 求平均后赋值
end
% 更新参数估计值
theta_hat = theta_hat - alpha * grad_J;
end
% 比较结果(可以根据需要进一步分析和可视化结果)
% disp(['使用u(t)=1时最终参数估计值: ', num2str(theta_hat(1)), ' ', num2str(theta_hat(2)), ' ', num2str(theta_hat(3))]);
disp(['使用u(t)=sin(0.2t)+sin(t)时最终参数估计值: ', num2str(theta_hat(1)), ' ', num2str(theta_hat(2)), ' ', num2str(theta_hat(3))]);