系统辨识副本

1、C+目录1、 LabWork12、 LabWork23、 LabWork34、 实验心得与体会5、 附录一（数据部分）6、 附录二（代码部分）01、 LabWork11. 请编出矩阵A与B相乘得到矩阵R的运算计算机程序。要求：(1) A和B的维数及数值可通过键盘及数据文件输入。 (2)计算结果R可由屏幕及文件输出。运行及结果显示部分-（程序见附录）文件中矩阵的存储格式：2 3 /（A矩阵维数2×3）1 2 3 3 4 53 2 /（B矩阵维数3×2）1 21 32 3-2. 矩阵求逆部分：运行及结果显示部分2、 LabWork2编写并调试动态模型仿真程序:模型：y(k)-1

2、.5y(k-1)+0.7y(k-2)=u(k-1)+0.5u(k-2)+v(k)已知白噪声v(k) 数据文件为DV, 数据长度L=500要求：(1) 产生长度为L的M 序列数据文件DU (2) 产生长度为L的模型输出数据文件DY运行及结果显示部分- 至此DU、DY数据已存储于硬盘中，可任意调用、读取。3、 LabWork3 & 4 编写并调试动态离散时间模型LS成批算法程序。要求：(1)原始数据由DU和DY读出； (2)调用求逆及相乘子程序； (3)显示参数辨识结果。-程序运行界面：- 编写并调试动态离散时间模型LS递推算法程序。要求：(1)原始数据由DU和DY读出； (2)显示辨识结

3、果； (3)设置选择变量决定是否输出中间结果。-程序运行界面：4、 实验心得与体会第一题矩阵乘法求解思路：对于矩阵相乘，我们并不陌生，主要是矩阵A对应行的值依次乘以矩阵B对应列的值，在最后把前面相乘的结果依次相加，其余的位置值的计算以此类推。我们最后就可以获得矩阵C的结果。根据线性代数定义中所说：当矩阵A的列数与矩阵B的行数不等时，矩阵乘法没有意义。故，我们要在矩阵输入后，判断矩阵的维数是否符合条件。矩阵结果从文件中读入和读出，我们主要用到C+的一个库函数：fstream 。这个函数比C语言中的FILE函数好用得多，它的文件的读写操作和iostream中的cout、cin操作极其相近。第一题矩

4、阵求逆求解思路：有根据公式：其中，|A|为行列式的值，为矩阵A的伴随矩阵。一开始，这部分的求解方法先考虑到了用代数余子式的方式求解。计算矩阵的逆，我们首先要做的就是求出矩阵的行列式（假设矩阵的行列式值不为零），然后再求出其伴随矩阵，最后两者按照公式相除，就得到了矩阵的逆。第二题求解思路：M序列是一种伪随机二进制序列。在实际中，我们需要产生周期足够长的序列，保证模拟现实干扰的情况。其实这个部分“最难”做的也就是M序列的产生。一开始本着利用课件里丰富的资源，如，用FORTRAN语言写的M序列产生函数，以及M序列的理论部分。想根据以上信息自行设计并产生相应的M序列数据。最后找了一种M序列的产生方法：

5、移位寄存器。这个方法很好的产生了一定周期的函数，而且也很好的符合了所需模拟的过程。在DY数据获取的部分，将其设置了初始值，Y(1)=0，Y(2)=0，初值的设定不对以后的辨识结果有太大干扰。最小二乘的结果足够精确的前提是数据的量远远大于未知数，即需要辨识的数据的量，就可以很好的辨识出系统的参数。最后将其余部分代入模型就可以依次算出Y的各个数据。第三题BLS、第四题RLS求解思路：在第三题，我们已经得到了辨识所需的算法BLS，再根据之前所写的矩阵求逆和矩阵乘法程序，只要在X、Y矩阵的生成上不出现问题，可以很轻松的完成这样一步工作。第四题，根据所给的RLS算法，根据以下的流程，我们可以很轻松的就计

6、算出所需辨识模型的参数。如果仅仅看问题部分，感觉就需要实现一个矩阵乘法、矩阵求逆、用移位寄存器原理实现M序列、由DV和DU得到DY、最后再套用两个最小二乘公式的问题。当然，在实际求解过程中，碰到了很多的问题。首先，一开始矩阵求逆的部分，就是因为原先线性代数没有学好的原因，错把三阶以下行列式值计算公式套用到了四阶以上。而实际上，四阶以上行列式的求解其实是要计算 n！，而不是原先错以为的8次乘法和加减法运算。其次，在获取DU、DY模型的时候，由于不是很清楚M序列的产生原理，一开始就不清楚要怎么做下去。最后找了一个MATLAB语言写的移位寄存器产生M序列的函数，并自己把它“翻译”成C+语言。之后在和

7、刘操他们交流后，发现用随机函数rand()也是一个可以产生M序列很好的方法。这个方法有一个好处也就是产生的序列周期比较长，而且内部的序列足够混乱，但一个不足就是在最后用最小二乘辨识参数的时候，辨识结果会比较随机。然后是X序列和Y序列的获取过程，这个部分因为在使用二维数组存储数据时未考虑好其空间维数，导致在在实际运行是发生了很多次数组越界的状况。解决方法是再次分析模型参数，解出两个未知量的值：n=4，m=498。其后的部分在此基础上求解就方便了很多。最后是算法的调用部分，一开始的BLS算法难度适中，主要是完成之前的乘法和求逆就可以进行计算。RLS相对来说就较复杂一些，需要分析出一开始k的初始值。

8、因为在C+中，数组的起始索引是从0开始，而实际模型中是从1开始，在编程时也需要进行一定的转换。还有一点就是在从X(k)中，k的解是3，但在实际分析时发现，这个X(3)的值其实就是BLS算法中用到的X矩阵的第一行数据，即k=0。在k的值确定后，调用递推公式进行计算，就比较快捷。在实际辨识过程中，发现辨识的结果与第二题的M序列值密切相关，即使周期已经加长到了。最后分析的结果如下：虽然周期足够长，但是计算其中的自相关系数，肯定与0不接近，由此就M序列，就对最终的辨识结果产生了极大的影响。这次的学习和编程实践是一件很快乐的事情，又增长了自己的见识，还协助他人解决了他们编程中的很多问题，同时自己也从他们

9、那里学习到了很多东西。主要是从大家解决问题的方案和思想中学到了很多，比如矩阵求逆部分还可以使用线性代数书里提及的初等行变换方法，移位寄存器产生M序列。五、附录一（部分数据）RLS算法的中间结果：六、附录二（代码部分）LabWork1#include "stdafx.h"#include <iostream>#include <vector>#include <fstream>using namespace std;typedef vector< vector<int>> Mat;/定义一个二维数组void Inpu

10、t(Mat &a, Mat &b); /数据的输入void Calculate(const Mat a,const Mat b, Mat &c); /矩阵运算void Output(const Mat c); /结果输出int _tmain(int argc, _TCHAR* argv)Mat a, b, c;Input(a, b);Calculate(a, b, c);Output(c);return 0;void Input(Mat &a, Mat &b)bool model;int row = 0, line = 0;cout<<&qu

11、ot;模式选择 ：0、从文件读取 1、手动输入"<<endl;cin>>model;if(model = 0)ifstream in("in.txt");in>> line>> row;int tmp = 0;for (int i = 0 ;i < line; +i)vector<int> temp1;for (int j = 0 ;j <row; +j)in>>tmp;temp1.push_back(tmp);a.push_back(temp1);in>> line

12、>> row;for (int i = 0 ;i < line; +i)vector<int> temp1;for (int j = 0 ;j <row; +j)in>>tmp;temp1.push_back(tmp);b.push_back(temp1);elsecout<<"请输入矩阵一行数和列数"<<endl;cin>>line >> row;int tmp = 0;for (int i = 0 ;i < line; +i)vector<int> temp

13、1;for (int j = 0 ;j <row; +j)cin>>tmp;temp1.push_back(tmp);a.push_back(temp1);cout<<"请输入矩阵二行数和列数"<<endl;cin>>line >> row;for (int i = 0 ;i < line; +i)vector<int> temp1;for (int j = 0 ;j <row; +j)cin>>tmp;temp1.push_back(tmp);b.push_back(te

14、mp1);void Calculate(const Mat a,const Mat b, Mat &c)int tmp = 0;for (int i = 0; i < a.size(); i+)vector<int> temp3;for (int j = 0; j < b0.size(); j+)for (int k = 0; k < b.size(); k+)tmp += aik * bkj;temp3.push_back(tmp);tmp = 0;c.push_back(temp3);void Output(const Mat c)cout<&l

15、t;"the answer is :"<<endl;for (int i = 0; i < c.size(); i+)for (int j = 0; j < c0.size(); j+)cout<< cij<<" "cout << endl;cout<<"是否将结果存之文件（y/n）"<<endl;ofstream out("Out.txt", ofstream:out);char str;cin>>str;if (st

16、r = 'y')for (int i = 0; i < c.size(); i+)for (int j = 0; j < c0.size(); j+)out<< cij<<" "out << "n"elsereturn;LabWork2#include "stdafx.h"#include "Data.h"using namespace std;int _tmain(int argc, _TCHAR* argv)vector<int> ve

17、c;int tmp;bool flag = true;while (flag)cout<<"请在此输入M序列的逻辑电路-9个由0、1组合的数值"<<endl;for (int i = 0; i < 9; i+)cin >> tmp;vec.push_back(tmp);for (int i = 0; i < vec.size(); +i)if (!(veci = 1 | veci = 0)flag = true;vec.clear(); /改动后新增内容cout << "请输入值为 0或1" &

18、lt;< endl;break;flag = false;CData MyData;MyData.DataTransform();MyData.M_Serial(vec);MyData.GetDY();cout << "Calculate Done !" << endl << "Open the Directory to get the result ." << endl;return 0;/stdafx.h的头文件（其他头文件声明全在此）#include <vector>#include

19、 <string>#include <iostream>#include <fstream>#include <iomanip>#include <cmath>using namespace std;/CData.h 的头文件#pragma once#include "stdafx.h"#define MAX_SIZE 100class CDatapublic:CData(void);CData(void);private:vector<double> m_DblData;vector<doubl

20、e> mseq;vector<double> m_DY_Data;public:void DataTransform(void);vector<double>& M_Serial(const vector<int> FbCon);void GetDY(void);/CData.cpp 的功能实现部分#include "StdAfx.h"#include "Data.h"CData:CData(void)mseq.push_back(0);m_DY_Data.push_back(0);m_DY_Data.p

21、ush_back(0);void CData:DataTransform(void)/将我们的输入文件(DV)转换一下，得到C+可以识别的数值ifstream in("DV");string tmp;while (! in.eof()in >> tmp;double digit = atof(tmp.c_str();m_DblData.push_back(digit);vector<double>& CData:M_Serial(const vector<int> FbCon)/注：这里的数组都是从1开始计数，下标为零0的没有参与

22、运算int n = FbCon.size();int N = pow(2.0, n) - 1;vector<double> register1(1,0);for (int i = 1; i <= n - 1; +i)register1.push_back(0);register1.push_back(1);mseq.push_back(1);for (int i = 2; i <= N; +i)vector<double> register2(2,0);int sum = 0;for (int k = 1; k <= n; +k)sum += regi

23、ster1.at(k) * FbCon.at(k-1);register2.at(1) = sum % 2;for (int j = 2; j <= n; +j)register2.push_back(register1j-1);/for (int j = 1; j <= n; +j)register1j = register2j;/mseq.push_back(register1n);ofstream out("DU.txt", ofstream:trunc);for (int i = 1; i < mseq.size() - 11; +i)out &l

24、t;< mseq.at(i) << "t"if (!(i % 5)out << endl;out.close();return mseq;void CData:GetDY(void)/从硬盘中读取已经存下的DU数据ifstream in("DU.txt");ofstream out("DY.txt", ofstream:trunc);double tmp;vector<double> My_DU;while (! in.eof()/数据输入in >> tmp;My_DU.push_

25、back(tmp);for (int i = 2; i < m_DblData.size(); +i)/公式计算部分double tmp_DY = 1.5 * m_DY_Datai-1 - 0.7 * m_DY_Datai-2 + My_DUi-1 + 0.5 * My_DUi-2 + m_DblDatai;m_DY_Data.push_back(tmp_DY);for (int i = 0; i < m_DblData.size(); +i)/数据输出到文件部分out <<setw(14) << setfill(' ') <<

26、m_DY_Datai ;if (!(i+1) % 5)out<<endl;out.close();LabWork3&4/stdafx.h的头文件（其他头文件声明全在此）#include <iostream>#include <fstream>#include <vector>#include <iomanip>#include <cmath>#include <algorithm>#include <string>#include <cstdlib>using namespace

27、 std;/#include "stdafx.h"#include "Matrix_.h"#include "Data.h"typedef vector<vector<double>> Mat; /定义一个二维向量void Get_DU_DY();/得到DU,DY的值void initial_(vector<double> &DY,vector<double> &DU, Mat &X_s);/初始化部分向量值void Use_BLS(vector<doubl

28、e> DY,vector<double> DU, Mat X_s);/BLS算法调用void Use_RLS(vector<double> DY,vector<double> DU, Mat X_s);/RLS算法调用Mat num2Mat(const double num); /数字到二维向量的转换CMatrix_ Tranposition(CMatrix_ matr); /矩阵转置void Median(const Mat matrix); /中间结果的输出显示部分int _tmain(int argc, _TCHAR* argv)Get_DU_D

29、Y();vector<double> DU,DY;Mat X_s;initial_(DY, DU, X_s);int num;cout << "请选择数字：1.BLS t2.RLS" <<endl;cin >> num;if(num = 1)Use_BLS(DY, DU, X_s);/BLS函数的调用：else if(num = 2)Use_RLS(DY, DU, X_s);/BLS函数的调用：elsecout << "Do Nothing." << endl;return 0;re

30、turn 0;void Get_DU_DY()/M序列部分代码/vector<int> vec;int tmpp = 0;cout << "推荐输入( 误差：<=1.5%) n"<< "010100101 n"<< "011110011" <<endl;cout << "或推荐输入( 误差：<=2.0%) n"<< "001011111 n"<< "010100101 n&quo

31、t;<< "011110011 n"<< "101100101 n"<< "101101001 n"<< "111000001 n"<< "111111001 n"<< "011110011" <<endl;cout << "输入逻辑电路值（0、1），共9位：" << endl;for (int k = 0; k < 9;k+)cin>&

32、gt;tmpp;vec.push_back(tmpp);CData MyData(0, 0);MyData.DataTransform();MyData.M_Serial(vec);MyData.GetDY();void initial_(vector<double> &DY,vector<double> &DU, Mat &X_s)/从文件中读取数据，上一次作业的ifstream in_DU("DU.txt");ifstream in_DY("DY.txt");ofstream out("sav

33、e.txt",ofstream:trunc);out.close();double tmp = 0;while(!in_DU.eof()in_DU>>tmp;DU.push_back(tmp);while(!in_DY.eof()in_DY>>tmp;DY.push_back(tmp);in_DU.close();in_DY.close();/X矩阵获取，对于BLS和RLS值都一样vector<double> tmp_val;int MAX_M_VALUE = 498;for (int i = 0; i < MAX_M_VALUE; i+)f

34、or (int j = i + 1; j >= i; j-)tmp_val.push_back(-1*DYj);for (int j = i + 1; j >= i; j-)tmp_val.push_back(DUj);X_s.push_back(tmp_val);tmp_val.clear();Mat num2Mat(const double num)Mat m;vector<double> tmp;tmp.push_back(num);m.push_back(tmp);return m;CMatrix_ Tranposition(CMatrix_ matr)/矩阵转

35、置CMatrix_ m;vector<double> tmp;for (newint i = 0; i < matr.m_Matrix0.size(); +i)for (newint j = 0; j < matr.m_Matrix.size(); +j)tmp.push_back(matr.m_Matrixji);m.m_Matrix.push_back(tmp);tmp.clear();return m;/以下是程序最核心的部分void Use_BLS(vector<double> DY, vector<double> DU, Mat X_s

36、) /BLS算法int MAX_M_VALUE = 498;/Y 矩阵的获得vector<double> get_Y;Mat Get_FiY;for (int i = 2; i < MAX_M_VALUE + 2; i+)get_Y.push_back(DYi);Get_FiY.push_back(get_Y);get_Y.clear();/CMatrix_ matrix_X(X_s);/类部分CMatrix_ matrix_Y(Get_FiY);/CMatrix_ LS_result;/调用BLSCMatrix_ tmp_aa;/tmp_aa = Tranposition(

37、matrix_X) * matrix_X;double tm_value = tmp_aa.det(tmp_aa.m_Matrix, tmp_aa.m_Matrix.size();tmp_aa = tmp_aa.Reverse(tmp_aa,tm_value);LS_result = tmp_aa * Tranposition(matrix_X) * matrix_Y;cout << "以下是运行结果：" << endl;LS_result.Print();cout<<endl;/void Use_RLS(vector<double

38、> DY,vector<double> DU, Mat X_s) /RLS算法int MAX_M_VALUE = 498;/P矩阵初始化Mat P_s;vector<double> tmp_val;for (int i = 0; i < 4; i+)for (int j = 0; j < 4; j+)if (i = j)tmp_val.push_back(1*1000000)elsetmp_val.push_back(0);P_s.push_back(tmp_val);tmp_val.clear();/thita 初始化Mat thita_s;for

39、(int j = 0; j < 4; j+)tmp_val.push_back(0);thita_s.push_back(tmp_val);tmp_val.clear();/存储中间值的类CMatrix_ matrix_Median;/CMatrix_ matrix_thita(thita_s);CMatrix_ matrix_P1(P_s);CMatrix_ matrix_gama;double gama = 0;/matrix_Median.m_Matrix.push_back(num2Mat(gama).at(0);matrix_Median.m_Matrix.push_back(

40、Tranposition(thita_s).m_Matrix.at(0);matrix_Median.m_Matrix.push_back(matrix_P1.m_Matrix.at(0);matrix_Median.m_Matrix.push_back(matrix_P1.m_Matrix.at(1);matrix_Median.m_Matrix.push_back(matrix_P1.m_Matrix.at(2);matrix_Median.m_Matrix.push_back(matrix_P1.m_Matrix.at(3);/for (int i = 1; i < MAX_M_V

41、ALUE; i+)CMatrix_ matrix_X_tmp(X_s.at(i-1);CMatrix_ matrix_X = Tranposition(matrix_X_tmp);CMatrix_ matrix_Yk(num2Mat(DY.at(i+1);/gamaCMatrix_ gama_tmp = Tranposition(matrix_X) * matrix_P1 * matrix_X;gama = (1.0/(gama_tmp.m_Matrix00 + 1);/thitamatrix_thita = matrix_thita + gama * matrix_P1 * matrix_X

42、 * (matrix_Yk - Tranposition(matrix_X) * matrix_thita);/Pmatrix_P1 = matrix_P1 - gama * matrix_P1 * matrix_X * Tranposition(matrix_X) * matrix_P1;/matrix_Median.m_Matrix.push_back(num2Mat(gama).at(0);matrix_Median.m_Matrix.push_back(Tranposition(matrix_thita).m_Matrix.at(0);matrix_Median.m_Matrix.pu

43、sh_back(matrix_P1.m_Matrix.at(0);matrix_Median.m_Matrix.push_back(matrix_P1.m_Matrix.at(1);matrix_Median.m_Matrix.push_back(matrix_P1.m_Matrix.at(2);matrix_Median.m_Matrix.push_back(matrix_P1.m_Matrix.at(3);/matrix_thita.Print();Median(matrix_Median.m_Matrix);void Median(const Mat matrix)int count =

44、 1;char m;cout<<"中间结果输出到'Median.txt'?（y/n）" << endl;cin >> m;if (m = 'y' | m = 'Y')ofstream out("Median.txt",ofstream:trunc);for (newint i = 0; i < matrix.size(); )out<<"第" << count+ << "次" <<

45、; "递推"<<endl;out<<"gama: " << matrixi0 <<endl;out<<"thita : "i+;for (newint j = 0; j < matrixi.size(); j+)out<<matrixij<<" "out<<"'" <<endl;i+;out<<"P 矩阵:" <<endl;for

46、(newint j = i; j < i + 4; j+)for (newint k = 0; k < matrixj.size(); k+)out <<setprecision(5)<<setw(13) << setfill(' ') << matrixjk;out<<endl;i = i + 4;out<<"/"<<endl;cout<<"已将结果输出到 'Median.txt'!"<<endl;ou

47、t.close();/CMatrix.h/#pragma oncetypedef vector<vector<double>> Mat;typedef vector<double>:size_type newint;class CMatrix_public:CMatrix_(void);CMatrix_(Mat matr);CMatrix_(vector<double> matr);CMatrix_(void);public:double det(Mat matrix, int size); /行列式值的计算（递归）CMatrix_ Revers

48、e(const CMatrix_ matrix, double value); /矩阵求逆/矩阵乘法的运算符重载friend CMatrix_ operator*(const CMatrix_ mat1, const CMatrix_ mat2);/乘法的运算符重载（第一位数字，第二位矩阵）friend CMatrix_ operator*(const double num, const CMatrix_ mat1);/矩阵加法的运算符重载friend CMatrix_ operator+(const CMatrix_ mat1, const CMatrix_ mat2);/矩阵减法的运算符重载friend CMatrix_ operator-(const CMatrix_ mat1, const CMatrix_ mat2);/矩阵间的赋值语句CMatrix_ operator=(const CMatrix_ mat2);/输出当前矩阵值void Print();public:Mat m_Mat