// // ARX 模型参数估计MATLAB estim / arx 的 C 实现// 编译 g -stdc17 -O2 arx_estim.cpp -o arx_estim// #includevector#includestring#includecmath#includenumeric#includealgorithm#includeiostream#includerandom// ------------------------------------------------------------// 数据结构// ------------------------------------------------------------structARXModel{intna;intnb;intnk;std::vectordoublea;// a1..anastd::vectordoubleb;// b1..bnbdoubleTs;ARXModel(intna_0,intnb_0,intnk_1,doubleTs_1.0):na(na_),nb(nb_),nk(nk_),a(na_,0.0),b(nb_,0.0),Ts(Ts_){}};structEstimResult{ARXModel model;doublelossFunction;doublefitPercent;boolsuccess;std::string message;};// ------------------------------------------------------------// 构建回归向量 phi(t) [ -y(t-1)..-y(t-na), u(t-nk)..u(t-nk-nb1) ]// ------------------------------------------------------------std::vectordoublebuildRegressor(conststd::vectordoubley,conststd::vectordoubleu,intt,intna,intnb,intnk){std::vectordoublephi(nanb,0.0);for(inti1;ina;i)phi[i-1]-y[t-i];for(inti0;inb;i)phi[nai]u[t-nk-i];returnphi;}// ------------------------------------------------------------// 求解正规方程 (Phi^T Phi) theta Phi^T Y// 使用高斯消元 部分主元返回残差平方和// ------------------------------------------------------------boolsolveNormalEquations(conststd::vectorstd::vectordoublePhi,conststd::vectordoubleY,std::vectordoubletheta,doubleresidualSS){intnstatic_castint(Phi[0].size());intNstatic_castint(Phi.size());std::vectorstd::vectordoubleA(n,std::vectordouble(n,0.0));std::vectordoublerhs(n,0.0);for(inti0;in;i){for(intj0;jn;j){doubles0.0;for(intk0;kN;k)sPhi[k][i]*Phi[k][j];A[i][j]s;}doublesy0.0;for(intk0;kN;k)syPhi[k][i]*Y[k];rhs[i]sy;}std::vectorintperm(n);std::iota(perm.begin(),perm.end(),0);for(intcol0;coln;col){intpivotcol;doublemaxAbsstd::abs(A[perm[col]][col]);for(intrcol1;rn;r){doublevstd::abs(A[perm[r]][col]);if(vmaxAbs){maxAbsv;pivotr;}}std::swap(perm[col],perm[pivot]);doublediagA[perm[col]][col];if(std::abs(diag)1e-12)returnfalse;for(intrcol1;rn;r){doublefA[perm[r]][col]/diag;for(intccol;cn;c)A[perm[r]][c]-f*A[perm[col]][c];rhs[perm[r]]-f*rhs[perm[col]];}}theta.assign(n,0.0);for(intin-1;i0;--i){doublesrhs[perm[i]];for(intji1;jn;j)s-A[perm[i]][j]*theta[j];theta[i]s/A[perm[i]][i];}residualSS0.0;for(intk0;kN;k){doublepred0.0;for(intj0;jn;j)predPhi[k][j]*theta[j];doubleeY[k]-pred;residualSSe*e;}returntrue;}// ------------------------------------------------------------// 主估计函数等价于 MATLAB 中的 arx / estim(...,arx)// ------------------------------------------------------------EstimResultestimateARX(conststd::vectordoubley,conststd::vectordoubleu,intna,intnb,intnk){EstimResult result;result.successfalse;intNstatic_castint(y.size());if(static_castint(u.size())!N){result.messagey 和 u 长度不一致;returnresult;}if(na0||nb0){result.messagena 和 nb 必须为正整数;returnresult;}intstartIdxstd::max(na,nknb-1);if(startIdxN){result.message数据长度不足;returnresult;}intMN-startIdx;intnParamsnanb;if(MnParams){result.message有效数据点少于参数个数;returnresult;}std::vectorstd::vectordoublePhi(M,std::vectordouble(nParams));std::vectordoubleY(M);for(intk0;kM;k){inttstartIdxk;Phi[k]buildRegressor(y,u,t,na,nb,nk);Y[k]y[t];}std::vectordoubletheta;doubleresidualSS0.0;if(!solveNormalEquations(Phi,Y,theta,residualSS)){result.message正规方程求解失败矩阵奇异;returnresult;}result.model.nana;result.model.nbnb;result.model.nknk;result.model.a.assign(theta.begin(),theta.begin()na);result.model.b.assign(theta.begin()na,theta.end());doublemeanY0.0;for(intkstartIdx;kN;k)meanYy[k];meanY/M;doubletotalSS0.0;for(intkstartIdx;kN;k){doubledy[k]-meanY;totalSSd*d;}result.lossFunctionresidualSS;result.fitPercent(totalSS1e-15)?(1.0-residualSS/totalSS)*100.0:0.0;result.successtrue;result.message估计成功;returnresult;}// ------------------------------------------------------------// 一步预测// ------------------------------------------------------------std::vectordoublepredictARX(constARXModelmodel,conststd::vectordoubley,conststd::vectordoubleu){intNstatic_castint(y.size());std::vectordoubleyPred(N,0.0);intstartIdxstd::max(model.na,model.nkmodel.nb-1);for(inttstartIdx;tN;t){doublepred0.0;for(inti1;imodel.na;i)pred-model.a[i-1]*y[t-i];for(inti0;imodel.nb;i)predmodel.b[i]*u[t-model.nk-i];yPred[t]pred;}returnyPred;}// ------------------------------------------------------------// 示例主程序// ------------------------------------------------------------intmain(){constintN500;std::vectordoubleu(N,0.0),y(N,0.0);std::mt19937rng(42);std::normal_distributiondoubledistU(0.0,1.0);std::normal_distributiondoubledistNoise(0.0,0.1);for(intt0;tN;t)u[t]distU(rng);// 真实系统: y(t) - 1.5y(t-1) 0.7y(t-2) 0.8u(t-1) 0.3u(t-2) e(t)// 即 a [-1.5, 0.7], b [0.8, 0.3], nk 1for(intt2;tN;t){y[t]1.5*y[t-1]-0.7*y[t-2]0.8*u[t-1]0.3*u[t-2]distNoise(rng);}autoresestimateARX(y,u,/*na*/2,/*nb*/2,/*nk*/1);if(res.success){std::cout 估计结果 \n;std::couta [;for(doublev:res.model.a)std::coutv ;std::cout]\n;std::coutb [;for(doublev:res.model.b)std::coutv ;std::cout]\n;std::cout拟合度: res.fitPercent%\n;std::cout损失函数: res.lossFunction\n;}else{std::cerr估计失败: res.message\n;}return0;}