控制系统程序变换的矩阵变换算法

!********************************************************
! 控制系统程序变换的矩阵变换算法
!********************************************************
! SINOMACH
!
! FORTRAN95,2013&18 version
!
! 控制系统程序变换（Program Transition）的
! 矩阵（Matrix Operation）变换算法
!
!
! 周勇, JohnBktJr, GoldbergJr,YidiW,ZhendongF,
! EricL, NadarJr, LilyZ, KnankJha
! 20221125
!********************************************************
!c $$$ Run, Output:
!********************************************************
!********************************************************
非时变控制系统设计和检测调试稳定性及可靠性中有参数和数据采集和分析
程序变换是一个数字计算机程序用于得到矩阵的变换，phi(t）的一
个线性时间不变量系统。程序还在用户选择的时间段内确定状态变量的
值。程序使用一个缩短的矩阵系列去确定 phi(t)=e**Ft*sigmak=0,N (Ft)**k/k!, (E.1), 这里
N+1 是使用的项的数目。当然， 当 N 等于无限时，这个指数系列准确地等于 phi(t). 当 N
是等于 10 或者 20 项时，一个缩短的系列常常可以满足常规的精度要求。可是，在每种情况下，
用户使用这个程序可以选择所要使用的项的个数 N。
这个程序考虑了线性向量微分方程
x'=Ax+Bu, (E.2)
这里 x=状态向量，及 u=输入信号。 一个新的向量 y 定义为
y=[x]
[m], (E.3)
这里 m 用于产生输入信号 u（t）， 并且状态向量和输入向量联接在一起作为一个向量。
那么这个新的向量微分方程是
y'=Fy. (E.4)
那么方程式（E.4）的解决方案是
y（t）=phi(t)y(0)
=eFty(0). (E.5)
为了估算在特定时间段的状态变量, 我们让 t=T,并且我们有
phi(t)=e**FT. (E.6)
缩短的系列使用 N+1 项, 在程序里标记为 IPROX 项, 所以
e**FT=~ I+ FT+ (FT)**2/2+ ...+ (FT)**IPROX/(IPROX)! (E.7)
必要渐进的项目数是精度要求的一个函数.可是要注意 1/10!=0.276*10**-6和1/20!=0.41*10-18.
通常15或者20项会足够,并且一个简单的方式监测精度是计算phi/(T),当 N=15 和 N=20 两者.
渐进对于一个小的项目数时准确,当选择了 T 时,以致 fijT<=1 给 F 的所有元素.
输入向量 m 产生所有多项式的输入信号或者正切(sine)波形输入.
所以一个阶跃或者斜坡输入会便于使用.输入向量用于产生系统输入 u(t),
所以对于一个单位斜坡输入,例如,我们有
m'=[m1'] [0 1] [m1]
= (E.8)
[m2'] [0 0] [m2]
并且
m(0)=[0,1]**T.
一般地我们使用次序 r 的输入向量 m=[m1,....mr]**T.
这个程序是以弗传(Fortran90,95,2003)语言编写. (n+r)可以有的最大值是
等于 13, 这里 n 是系统的次数(状态变量的数目).
使用了三组数据卡.首先列出了时间段 T, 值(n+r)=LIM, 以及缩
减系列 IPROX 的项目数.这个列表应该按照格式 F5.3,2I2. 第二组数
据的列表是 F 矩阵的元素. 元素应该按照格式 5E14.7 列出,以如下的
次序: f11, f12,...,f1LIM,f21,f22,...,f2LIM,...,fLIM1,fLIM2,...
,fLIMLIM.数据组列出主题状态向量的初始值,y(0).这些是根据格式
5E14.7 列表，次序为 x1(0),x2(0,...,xn(0),m1(0),m2(0),...,mr(0).
程序的输出。程序显示在这个附录上.
!*****************************************************************

R(s)--> O-->[5/s*(s+2)(s+4)]--*-- X1(s)

• ^| |
  |_____________________|
  图 E.1 （a）,

5      1/s  x3  1/s x2  1/s x1
u(t)=r(t) O-->---O--->---O--->---0--->---0          
                 |       |       |       |
                 |--<----|       |       |
                 |  -3           |       |
                 |------<--------|       |
                 |      -2               |
                 |----------<------------|
                            -5
         图 E.1（b）

例 E.1. 让我们考虑一下如图 E.1(a)所示的系统，这里状态变量辨
识如图 E.1(b).状态向量的微分方程是
[ 0 1 0] [0]
x'=[ 0 0 1]x+[0]r(t). (E.9)
[-5 -2 -3] [5]
考虑这个系统有一个单位的阶跃输入，这样 m1（t）=r（t），m1'.=0,及 m1（0）=1。
那么主题系统的矩阵方程式就是
y'=d/dt[x1/x2/x3/m1]=[0 1 0 0/0 0 1 0/-5 -2 -3 5/0 0 0 0][x1/x2/x3/m1] (E.10)
以及 y**T（0）=[0,0,0,1].
初始让我们取 T=0.5 和 N+1=IPROX=10 项。另外注意 LIM=（n+r）=4。
那么第一张数据清单是
0.5 0 0 0 0 4 0 1 0
|-----| |---| |---|
T LIM IPROX
F 矩阵的元素是：
0.5 0 0 0 0 4 0 1 0
0.0 E+00 1.0 E+00 0.0 E+00 0.0 E+00 0.0 E+00
0.0 E+00 1.0 E+00 0.0 E+00 -5.0 E+00 -2.0 E+00
-3.0 E+00 5.0 E+00 0.0 E+00 0.0 E+00 0.0 E+00
0.0 E+00
对 y（t）的初始条件是：
0.0 E+00 0.0 E+00 0.0 E+00 1.0 E+00
计算的phi(T)值和 t=kT 的状态变量如下：
PHI MATRIX ELEMENTS PHI(I,J) IN THE ORDER
PHI(1,1), ......................,PHI(LIM,LIM).
###############################
###############################
STATE VECTOR ELEMENTS IN THE ORDER
X(1),..................,X(N), M(1),.......,M(R).
###############################
###############################
................
当 IPROX 增加到 15，得到下列输出：
###############################
................
当 IPROX 增加到 20 时，得到下列 PHI 矩阵：
###############################
................
明显地， 在这种情况下 10 项就足够了， 15 项会提供五位数的精
度。 很清楚程序容易使用，项的数目可以方便地更换。
!********************************************************************
!********************************************************************
Program transition
Dimension a(13, 13), b(13, 13), c(13, 13), d(13, 13), x(13), y(13)
77 Read (*, 202) t, lim, iprox
Read (*, 201) ((a(i,j),j=1,lim), i=1, lim)
Read (*,201) (y(i), i=1, lim)
Do i = 1, lim
Do j = 1, lim
a(i, j) = t*a(i, j)
End Do
End Do
Do i = 1, lim
Do j = 1, lim
c(i, j) = 0.0
d(i, j) = a(i, j)
b(i, j) = a(i, j)
End Do
End Do
Do ik = 2, iprox
den = ik
Do i = 1, lim
Do j = 1, lim
Do k = 1, lim
c(i, j) = c(i, j) + a(i, k)*b(k, j)/den
End Do
End Do
End Do
Do i = 1, lim
Do j = 1, lim
b(i, j) = c(i, j)
d(i, j) = d(i, j) + c(i, j)
End Do
End Do
Do i = 1, lim
Do j = 1, lim
c(i, j) = 0.0
End Do
End Do
End Do
Do i = 1, lim
c(i, i) = 1.0
End Do
Do i = 1, lim
Do j = 1, lim
d(i, j) = d(i, j) + c(i, j)
End Do
End Do
Write (*, 29)
29 format('PHI MATRIX ELEMENTS PHI(I,J) IN THE ORDER', /,&
'phi(1,1),phi(1,2),phi(1,3),phi(1,4),phi(1,5),phi(1,6),phi(1,7),&
phi(1,8), phi(1,9), phi(1,10), phi(1,11), phi(1,12), phi(1,lim),&
phi(2,1), phi(2,2), phi(2,3), phi(2,4), phi(2,5), phi(2,6), phi(2,7),&
phi(2,8), phi(2,9), phi(2,10), phi(2,11), phi(2,12), phi(2,lim), /,&
phi(3,1), phi(3,2), phi(3,3), phi(3,4), phi(3,5), phi(3,6), phi(3,7),&
phi(3,8), phi(3,9), phi(3,10), phi(3,11), phi(3,12), phi(3,lim),/,&
phi(4,1), phi(4,2), phi(4,3), phi(4,4), phi(4,5), phi(4,6), phi(4,7),&
phi(4,8), phi(4,9), phi(4,10), phi(4,11), phi(4,12), phi(4,lim),/,&
phi(5,1), phi(5,2), phi(5,3), phi(5,4), phi(5,5), phi(5,6), phi(5,7),&
phi(5,8), phi(5,9), phi(5,10), phi(5,11), phi(5,12), phi(5,lim),/,&
phi(6,1), phi(6,2), phi(6,3), phi(6,4), phi(6,5), phi(6,6), phi(6,7),&
phi(6,8), phi(6,9), phi(6,10), phi(6,11), phi(6,12), phi(6,lim),/,&
phi(7,1), phi(7,2), phi(7,3), phi(7,4), phi(7,5), phi(7,6), phi(7,7),&
phi(7,8), phi(7,9), phi(7,10), phi(7,11), phi(7,12), phi(7,lim),/,&
phi(8,1), phi(8,2), phi(8,3), phi(8,4), phi(8,5), phi(8,6), phi(8,7),&
phi(8,8), phi(8,9), phi(8,10), phi(8,11), phi(8,12), phi(8,lim),/,&
phi(9,1), phi(9,2), phi(9,3), phi(9,4), phi(9,5), phi(9,6), phi(9,7),&
phi(9,8), phi(9,9), phi(9,10), phi(9,11), phi(9,12), phi(9,lim),/,&
phi(10,1), phi(10,2), phi(10,3), phi(10,4), phi(10,5), phi(10,6),phi(10,7),&
phi(10,8), phi(10,9), phi(10,10), phi(10,11), phi(10,12),phi(10,lim),/,&
phi(11,1), phi(11,2), phi(11,3), phi(11,4), phi(11,5),phi(11,6), phi(11,7),&
phi(11,8), phi(11,9), phi(11,10), phi(11,11),phi(11,12), phi(11,lim),/,&
phi(12,1), phi(12,2), phi(12,3), phi(12,4),phi(12,5), phi(12,6), phi(12,7),&
phi(12,8), phi(12,9), phi(12,10),phi(12,11), phi(12,12), phi(12,lim),/,&
phi(lim,1), phi(lim,2),phi(lim,3), phi(lim,4), phi(lim,5), phi(lim,6), &
phi(lim,7),phi(lim,8), phi(lim,9), phi(lim,10), phi(lim,11),&
phi(lim,12),phi(lim,lim).')
Write(*,30)((d(i,j),j=1,lim),i=1,lim)
Write(*,203)
203 format('STATE VECTOR ELEMENTS IN THE ORDER',/,&
'x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),x(9),x(10),x(11),&
x(12),x(n),m(1),m(2),m(3),m(4),m(5),m(6),m(7),m(8),m(9),m(10),&
m(11),m(12),m(r) ')
Do l=1,00
Do i=1,lim
x(i)=0.0
End Do
Do i=1,lim
Do j=1,lim
x(i)=x(i)+d(i,j)*y(j)
End Do
End Do
Do k=1,lim
y(k)=x(k)
End Do
Write(*,30)(y(i),i=1,lim)
End Do
Goto 77
202 Format(F5.3,2I3)
201 Format(5E15.7)
30 Format(4(4X,E15.9))
End Program transition
!**********************************
!**********************************
!c $$
Standard Input:
!**********************************
!**********************************
0.5,4,10
0.0,1.0,0.0,0.0,0.0,
0.0,1.0,0.0,-5.0,-2.0,
-3.0,5.0,0.0,0.0,0.0,
0.0
0.0,0.0,0.0,1.0
0.5, 4, 10
0.0,1.0,0.0,0.0,0.0,
0.0,1.0,0.0,-5.0,-2.0,
-3.0,5.0,0.0,0.0,0.0,
0.0
0.2,0.3,0.4,1.0
0.6, 4, 5
0.0,1.0,0.0,0.0,0.0
0.0,1.0,0.0,-5.0,-2.0
-3.0,5.0,0.0,0.0,0.0
0.0
0.0,6.0,0.0,2.0
0.7, 4, 6
0.0,1.0,0.0,0.0,0.0
0.0,1.0,0.0,-5.0,-2.0
-3.0,5.0,0.0,0.0,0.0
0.0
0.0,0.0,0.0,3.0
0.8, 4, 12
0.0,1.0,0.0,0.0,0.0
0.0,1.0,0.0,-5.0,-2.0
-3.0,5.0,0.0,0.0,0.0
0.0
0.0,0.0,4.0,4.0
0.9, 4, 15
0.0,1.0,0.0,0.0,0.0
0.0,1.0,0.0,-5.0,-2.0
-3.0,5.0,0.0,0.0,0.0
0.0
0.0,3.0,0.0,5.0
1.0, 4, 20
0.0,1.0,0.0,0.0,0.0
0.0,1.0,0.0,-5.0,-2.0
-3.0,5.0,0.0,0.0,0.0
0.0
2.0,0.0,0.0,6.0
0.1, 4, 13
0.0,1.0,0.0,0.0,0.0
0.0,1.0,0.0,-5.0,-2.0
-3.0,5.0,0.0,0.0,0.0
0.0
1.0,0.5,0.0,1.0
!********************************************************
!********************************************************
!c