public_ident.bas

系统辨识

Attribute VB_Name = "Public_Ident"
' 本工程功能为：递推最小二乘法
'待辨识的 模型2 结构为：
'y(t)+a1y(t-1)+a2y(t-2)=b1u(t-1)+b2u(t-2)
'该程序参考《控制系统的计算机辅助设计程序汇编》
'清华出版 孙增圻(qi),并略有改变
'改变后的特点：速度比原来快一倍,所需变量比原来少22个（模型为2阶时）
'由文件 Public_Ident.bas 和 frmMain.frm组成
' 下面是文件Public_Ident.bas的源代码
'参照公式:《自适应控制系统》P73 谢新民 清华出版
Option Explicit
'[注释]  辨识所需要的公共变量
Public Const rou = 0.995 '遗忘因子
Public Const Order = 2 '模型的阶次
Public Const P_Len = Order + Order '矩阵P的维数
Public Seita_t(P_Len) As Double  't时刻辨识的参数
Public Model_Error As Double         '记录辨识的模型的残差平方和

Sub Ident(Ini_flag As Boolean, Control_Ut As Double, Samp_Yt As Double)
'注释：Ini_flag = True '表示调用该函数进行辨识    初始化
'注释：Ini_flag = False '表示调用该函数进行       递推辨识
'注释：Control_Ut为 t 时刻的控制输入Um(t)
'注释：Samp_Yt为    t 时刻的系统输出Ym(t)

Static Ym(Order) As Double      '必须定义为静态变量或全局变量
Static Um(Order) As Double
Static Pt(P_Len, P_Len) As Double

Dim i, j, k  As Integer  '记录 FOR 循环次数
Dim Temp As Double
Dim Seita_t1(P_Len) As Double  't+1时刻辨识的参数
Dim Xt1(P_Len) As Double                '存放X(t+1)=[-y(t),-y(t-1),u(t+1),u(t),u(t-1)]
Dim epsilong As Double   '预测偏差epsilong = Samp_Yt -   X(t+1) * Seita_t

Dim Lt1(P_Len) As Double
Dim Lt1_Xt1(P_Len, P_Len) As Double
Dim Lt1_Xt1_Pt(P_Len, P_Len) As Double


If Ini_flag = True Then
'******************** 初始化开始 **********************
    '初始化Seita_t = [ -y(t-1), -y(t-2), u(t), u(t-1), u(t-2) ]
    For i = 1 To P_Len   '设参数的初始值均为0
    Seita_t(i) = 0.0001
    Next i
    '初始化矩阵  P
    For i = 1 To P_Len
    For j = 1 To P_Len
      Pt(i, j) = 0
    Next j
    Pt(i, i) = 10000000000#
    Next i
    '初始化X(t)
    '将X(t)向前移动一步,变成X(t+1)
    For i = Order To 1 Step -1
    Ym(i) = Ym(i - 1)
    Um(i) = Um(i - 1)
    Next i

    Ym(0) = Samp_Yt
    Um(0) = Control_Ut
   '初始化模型残差平方和
   Model_Error = 0
'******************** 初始化结束 **********************
Else


'********************  辨识 开始 **********************

'将X(t)向前移动一步,变成X(t+1)
For i = Order To 1 Step -1
Ym(i) = Ym(i - 1)
Um(i) = Um(i - 1)
Next i

Ym(0) = Samp_Yt
Um(0) = Control_Ut

For i = 1 To Order
Xt1(i) = -Ym(i)
Next i

For i = Order + 1 To P_Len
Xt1(i) = Um(i - Order)
Next i



'注释：Debug.Print观察送入的数据是否正确，对程序无影响，可删除
'Debug.Print Xt1(1); Xt1(2); Xt1(3); Xt1(4); Xt1(5)

'计算L(t+1)
For i = 1 To P_Len
Lt1(i) = 0
For j = 1 To P_Len
Lt1(i) = Lt1(i) + Pt(i, j) * Xt1(j)
Next j
Next i

Temp = 0
For i = 1 To P_Len
Temp = Temp + Xt1(i) * Lt1(i)
Next i
Temp = Temp + rou

For i = 1 To P_Len
Lt1(i) = Lt1(i) / Temp
Next i

'计算Seita(t+1)
Temp = 0
For j = 1 To P_Len
  Temp = Temp + Xt1(j) * Seita_t(j)
Next j
epsilong = Samp_Yt - Temp
  
For i = 1 To P_Len
  Seita_t1(i) = Seita_t(i) + Lt1(i) * epsilong
Next i

For i = 1 To P_Len
Seita_t(i) = Seita_t1(i)
Next i


'计算P(t+1)
For i = 1 To P_Len
For j = 1 To P_Len
Lt1_Xt1(i, j) = Lt1(i) * Xt1(j)
Next j
Next i


For i = 1 To P_Len
For j = 1 To P_Len
Temp = 0
For k = 1 To P_Len
Temp = Temp + Lt1_Xt1(i, k) * Pt(k, j)
Next k
Lt1_Xt1_Pt(i, j) = Temp
Next j
Next i

For i = 1 To P_Len
   For j = 1 To P_Len
     Pt(i, j) = (Pt(i, j) - Lt1_Xt1_Pt(i, j)) / rou
   Next j
Next i

'计算残差
Model_Error = Model_Error * rou + epsilong * epsilong
'********************  辨识 结束 **********************
End If

End Sub