丢番图（Diophantine）方程MATLAB求解
丢番图 (Diophantine) 方程在多项式中的一般形式为
A
(
z
−
1
)
X
(
z
−
1
)
+
B
(
z
−
1
)
Y
(
z
−
1
)
=
C
(
z
−
1
)
A(z^{-1}) X(z^{-1})+B(z^{-1}) Y(z^{-1})=C(z^{-1})
A(z−1)X(z−1)+B(z−1)Y(z−1)=C(z−1)
式中,
A
(
z
−
1
)
,
B
(
z
−
1
)
A(z^{-1}), B\left(z^{-1}\right)
A(z−1),B(z−1) 和
C
(
z
−
1
)
C\left(z^{-1}\right)
C(z−1) 为三个已知非零多项式, 且
A
(
z
−
1
)
,
B
(
z
−
1
)
A\left(z^{-1}\right), B\left(z^{-1}\right)
A(z−1),B(z−1) 和
C
(
z
−
1
)
C\left(z^{-1}\right)
C(z−1) 三者之 间无公因子,
X
(
z
−
1
)
X\left(z^{-1}\right)
X(z−1) 和
Y
(
z
−
1
)
Y\left(z^{-1}\right)
Y(z−1) 为末知的待求多项式。
这是一个末知数个数多于方程个数的方程, 因此是一个不定方程。如果
A
(
z
−
1
)
A\left(z^{-1}\right)
A(z−1) 和
B
(
z
−
1
)
B\left(z^{-1}\right)
B(z−1) 互 质，也就是
A
(
z
−
1
)
A\left(z^{-1}\right)
A(z−1) 和
B
(
z
−
1
)
B\left(z^{-1}\right)
B(z−1) 之间没有公因子, 则上述方程总存在解, 并且存在多个解。例如, 如果
X
0
(
z
−
1
)
X_{0}\left(z^{-1}\right)
X0(z−1) 和
Y
0
(
z
−
1
)
Y_{0}\left(z^{-1}\right)
Y0(z−1) 是上述方程的一个特解, 则
X
(
z
−
1
)
=
X
0
(
z
−
1
)
−
B
(
z
−
1
)
Q
(
z
−
1
)
Y
(
z
−
1
)
=
Y
0
(
z
−
1
)
+
A
(
z
−
1
)
Q
(
z
−
1
)
\begin{aligned} &X\left(z^{-1}\right)=X_{0}\left(z^{-1}\right)-B\left(z^{-1}\right) Q\left(z^{-1}\right) \\ &Y\left(z^{-1}\right)=Y_{0}\left(z^{-1}\right)+A\left(z^{-1}\right) Q\left(z^{-1}\right) \end{aligned}
X(z−1)=X0(z−1)−B(z−1)Q(z−1)Y(z−1)=Y0(z−1)+A(z−1)Q(z−1)
也是该方程的解, 其中
Q
(
z
−
1
)
Q\left(z^{-1}\right)
Q(z−1) 为任意多项式。控制问题中我们通常是求出
X
(
z
−
1
)
X\left(z^{-1}\right)
X(z−1) 或
Y
(
z
−
1
)
Y\left(z^{-1}\right)
Y(z−1) 的 最小阶解, 而最小阶解则是 Diophantine 方程在最小阶限制下的唯一解。以确定
X
(
z
−
1
)
X\left(z^{-1}\right)
X(z−1) 的最 小阶解为例加以说明。
确定
X
(
z
−
1
)
X\left(z^{-1}\right)
X(z−1) 和
Y
(
z
−
1
)
Y\left(z^{-1}\right)
Y(z−1) 的阶次
n
X
n_{X}
nX 和
n
Y
n_{Y}
nY 的方法是根据方程两边
z
−
1
z^{-1}
z−1 的同次幂系数相等的原 则建立一组线性方程。因为
X
(
z
−
1
)
X\left(z^{-1}\right)
X(z−1) 和
Y
(
z
−
1
)
Y\left(z^{-1}\right)
Y(z−1) 的系数末知待求, 所以方程组中的末知数的个数 为
X
(
z
−
1
)
X\left(z^{-1}\right)
X(z−1) 和
Y
(
z
−
1
)
Y\left(z^{-1}\right)
Y(z−1) 的阶次
n
X
n_{X}
nX 和
n
Y
n_{Y}
nY 之和再加 2 , 即末知数个数为
n
X
+
n
Y
+
2
n_{X}+n_{Y}+2
nX+nY+2; 而方程组中方程 式的个数为
n
C
+
1
n_{C}+1
nC+1, 或
n
X
+
n
A
+
1
n_{X}+n_{A}+1
nX+nA+1, 或
n
Y
+
n
B
+
1
n_{Y}+n_{B}+1
nY+nB+1, 其中
n
A
,
n
B
n_{A}, n_{B}
nA,nB 和
n
C
n_{C}
nC 分别代表
A
(
z
−
1
)
,
B
(
z
−
1
)
A\left(z^{-1}\right), B\left(z^{-1}\right)
A(z−1),B(z−1) 和
C
(
z
−
1
)
C\left(z^{-1}\right)
C(z−1) 的阶次。我们可以看出丢番图方程有唯一解的条件是方程个数等于末知数个数, 即
n
X
+
n
Y
+
2
=
n
Y
+
n
B
+
1
n
X
+
n
Y
+
2
=
n
X
+
n
A
+
1
n
X
+
n
Y
+
2
=
n
C
+
1
\begin{aligned} &n_{X}+n_{Y}+2=n_{Y}+n_{B}+1 \\ &n_{X}+n_{Y}+2=n_{X}+n_{A}+1 \\ &n_{X}+n_{Y}+2=n_{C}+1 \end{aligned}
nX+nY+2=nY+nB+1nX+nY+2=nX+nA+1nX+nY+2=nC+1
确定
X
(
z
−
1
)
X\left(z^{-1}\right)
X(z−1) 的最小阶解, 也就是限制
n
X
<
n
B
n_{X}<n_{B}
nX<nB 的解, 在此条件下, 可以确定
n
X
=
n
B
−
1
n
Y
=
max
{
n
A
−
1
,
n
C
−
n
B
}
\begin{gathered} n_{X}=n_{B}-1 \\ n_{Y}=\max \left\{n_{A}-1, n_{C}-n_{B}\right\} \end{gathered}
nX=nB−1nY=max{nA−1,nC−nB}
若限制
n
C
≤
n
A
+
n
B
−
1
n_{C} \leq n_{A}+n_{B}-1
nC≤nA+nB−1, 则
n
Y
=
n
A
−
1
n_{Y}=n_{A}-1
nY=nA−1
以单容水箱控制器设计模型为例进行丢番图方程的求解：
A
(
z
−
1
)
y
(
k
)
=
z
−
d
B
(
z
−
1
)
u
(
k
)
A\left(z^{-1}\right) y(k)=z^{-d} B\left(z^{-1}\right) u(k)
A(z−1)y(k)=z−dB(z−1)u(k)
其中
A
(
z
−
1
)
=
1
−
0.7652
z
−
1
−
0.2297
z
−
2
B
(
z
−
1
)
=
−
0.0006
+
0.0048
z
−
1
d
=
1
\begin{gathered} A\left(z^{-1}\right)=1-0.7652 z^{-1}-0.2297 z^{-2} \\ B\left(z^{-1}\right)=-0.0006+0.0048 z^{-1} \\ d=1 \end{gathered}
A(z−1)=1−0.7652z−1−0.2297z−2B(z−1)=−0.0006+0.0048z−1d=1
选取闭环极点位置为
{
0.6225
+
0.1380
i
,
0.6225
−
0.1380
i
}
\{0.6225+0.1380 i, 0.6225-0.1380 i\}
{0.6225+0.1380i,0.6225−0.1380i}
则理想闭环极点方程为
T
(
z
−
1
)
=
1
−
1.2451
z
−
1
+
0.4066
z
−
2
,
T\left(z^{-1}\right)=1-1.2451 z^{-1}+0.4066 z^{-2},
T(z−1)=1−1.2451z−1+0.4066z−2,
采用极点配置控制方法来实现对单容水箱液位的控制, 对应的控制律为
H
(
z
−
1
)
u
(
k
)
=
E
(
z
−
1
)
w
(
k
)
−
G
(
z
−
1
)
y
(
k
)
H\left(z^{-1}\right) u(k)=E\left(z^{-1}\right) w(k)-G\left(z^{-1}\right) y(k)
H(z−1)u(k)=E(z−1)w(k)−G(z−1)y(k)
其中
H
(
z
−
1
)
、
E
(
z
−
1
)
H\left(z^{-1}\right) 、 E\left(z^{-1}\right)
H(z−1)、E(z−1) 和
G
(
z
−
1
)
G\left(z^{-1}\right)
G(z−1) 为多项式,
w
(
k
)
w(k)
w(k) 为参考输入。
采用保留过程所有零点的极点配置方法有：
采用保留过程全部零点的极点配置算法, 此时的极点配置方程应写成
A
(
z
−
1
)
H
(
z
−
1
)
+
B
(
z
−
1
)
G
(
z
−
1
)
=
T
(
z
−
1
)
A\left(z^{-1}\right) H\left(z^{-1}\right)+B\left(z^{-1}\right) G\left(z^{-1}\right)=T\left(z^{-1}\right)
A(z−1)H(z−1)+B(z−1)G(z−1)=T(z−1)
其中
B
(
z
−
1
)
=
z
−
d
B
(
z
−
1
)
B\left(z^{-1}\right)=z^{-d} B\left(z^{-1}\right)
B(z−1)=z−dB(z−1)
阶次限制关系为
n
H
=
n
B
−
1
n
G
=
n
A
−
1
n
T
≤
n
A
+
n
B
−
1
\begin{aligned} &n_{H}=n_{B}-1 \\ &n_{G}=n_{A}-1 \\ &n_{T} \leq n_{A}+n_{B}-1 \end{aligned}
nH=nB−1nG=nA−1nT≤nA+nB−1
由于
A
(
z
−
1
)
A\left(z^{-1}\right)
A(z−1) 与
B
(
z
−
1
)
B\left(z^{-1}\right)
B(z−1) 互质, 上式一定有解。
根据多项式系数进行编程：
syms zn g0 g1 h0 h1;
A = 1 - 0.7652*zn -0.2297*zn^2;
B = -0.0006*zn+0.0048*zn^2;
T = 1 - 1.2451*zn +0.4066*zn^2;
G = g0-g1*zn; % ng = nb - 1 = 1阶
H = h0-h1*zn; % nf = na - 1 = 1阶
% Calculate polynomial
poly = A*H+B*G-T;
[g0v, g1v, h0v,h1v] = solve(coeffs(poly, zn));
double([g0v, g1v, h0v,h1v])
当不加积分器时, 取
E
(
1
)
=
T
(
1
)
/
B
(
1
)
E(1)=T(1) / B(1)
E(1)=T(1)/B(1)
得到:
E
(
z
−
1
)
=
38.7119
,
G
(
z
−
1
)
=
59.5098
−
21.4817
z
−
1
,
H
(
z
−
1
)
=
1
−
0.445
z
−
1
,
\begin{gathered} E\left(z^{-1}\right)=38.7119, \\ G\left(z^{-1}\right)=59.5098-21.4817 z^{-1}, \\ H\left(z^{-1}\right)=1-0.445 z^{-1}, \end{gathered}
E(z−1)=38.7119,G(z−1)=59.5098−21.4817z−1,H(z−1)=1−0.445z−1,
2925
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
