企业官网建设流程全解析

// ============================================================// ARX 模型参数估计（MATLAB estim / arx 的 C++ 实现）// 编译： g++ -std=c++17 -O2 arx_estim.cpp -o arx_estim// ============================================================#include<vector>#include<string>#include<cmath>#include<numeric>#include<algorithm>#include<iostream>#include<random>// ------------------------------------------------------------// 数据结构// ------------------------------------------------------------structARXModel{intna;intnb;intnk;std::vector<double>a;// a1..anastd::vector<double>b;// 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-nb+1) ]// ------------------------------------------------------------std::vector<double>buildRegressor(conststd::vector<double>&y,conststd::vector<double>&u,intt,intna,intnb,intnk){std::vector<double>phi(na+nb,0.0);for(inti=1;i<=na;++i)phi[i-1]=-y[t-i];for(inti=0;i<nb;++i)phi[na+i]=u[t-nk-i];returnphi;}// ------------------------------------------------------------// 求解正规方程 (Phi^T Phi) theta = Phi^T Y// 使用高斯消元 + 部分主元，返回残差平方和// ------------------------------------------------------------boolsolveNormalEquations(conststd::vector<std::vector<double>>&Phi,conststd::vector<double>&Y,std::vector<double>&theta,double&residualSS){intn=static_cast<int>(Phi[0].size());intN=static_cast<int>(Phi.size());std::vector<std::vector<double>>A(n,std::vector<double>(n,0.0));std::vector<double>rhs(n,0.0);for(inti=0;i<n;++i){for(intj=0;j<n;++j){doubles=0.0;for(intk=0;k<N;++k)s+=Phi[k][i]*Phi[k][j];A[i][j]=s;}doublesy=0.0;for(intk=0;k<N;++k)sy+=Phi[k][i]*Y[k];rhs[i]=sy;}std::vector<int>perm(n);std::iota(perm.begin(),perm.end(),0);for(intcol=0;col<n;++col){intpivot=col;doublemaxAbs=std::abs(A[perm[col]][col]);for(intr=col+1;r<n;++r){doublev=std::abs(A[perm[r]][col]);if(v>maxAbs){maxAbs=v;pivot=r;}}std::swap(perm[col],perm[pivot]);doublediag=A[perm[col]][col];if(std::abs(diag)<1e-12)returnfalse;for(intr=col+1;r<n;++r){doublef=A[perm[r]][col]/diag;for(intc=col;c<n;++c)A[perm[r]][c]-=f*A[perm[col]][c];rhs[perm[r]]-=f*rhs[perm[col]];}}theta.assign(n,0.0);for(inti=n-1;i>=0;--i){doubles=rhs[perm[i]];for(intj=i+1;j<n;++j)s-=A[perm[i]][j]*theta[j];theta[i]=s/A[perm[i]][i];}residualSS=0.0;for(intk=0;k<N;++k){doublepred=0.0;for(intj=0;j<n;++j)pred+=Phi[k][j]*theta[j];doublee=Y[k]-pred;residualSS+=e*e;}returntrue;}// ------------------------------------------------------------// 主估计函数：等价于 MATLAB 中的 arx / estim(...,'arx')// ------------------------------------------------------------EstimResultestimateARX(conststd::vector<double>&y,conststd::vector<double>&u,intna,intnb,intnk){EstimResult result;result.success=false;intN=static_cast<int>(y.size());if(static_cast<int>(u.size())!=N){result.message="y 和 u 长度不一致";returnresult;}if(na<=0||nb<=0){result.message="na 和 nb 必须为正整数";returnresult;}intstartIdx=std::max(na,nk+nb-1);if(startIdx>=N){result.message="数据长度不足";returnresult;}intM=N-startIdx;intnParams=na+nb;if(M<nParams){result.message="有效数据点少于参数个数";returnresult;}std::vector<std::vector<double>>Phi(M,std::vector<double>(nParams));std::vector<double>Y(M);for(intk=0;k<M;++k){intt=startIdx+k;Phi[k]=buildRegressor(y,u,t,na,nb,nk);Y[k]=y[t];}std::vector<double>theta;doubleresidualSS=0.0;if(!solveNormalEquations(Phi,Y,theta,residualSS)){result.message="正规方程求解失败（矩阵奇异）";returnresult;}result.model.na=na;result.model.nb=nb;result.model.nk=nk;result.model.a.assign(theta.begin(),theta.begin()+na);result.model.b.assign(theta.begin()+na,theta.end());doublemeanY=0.0;for(intk=startIdx;k<N;++k)meanY+=y[k];meanY/=M;doubletotalSS=0.0;for(intk=startIdx;k<N;++k){doubled=y[k]-meanY;totalSS+=d*d;}result.lossFunction=residualSS;result.fitPercent=(totalSS>1e-15)?(1.0-residualSS/totalSS)*100.0:0.0;result.success=true;result.message="估计成功";returnresult;}// ------------------------------------------------------------// 一步预测// ------------------------------------------------------------std::vector<double>predictARX(constARXModel&model,conststd::vector<double>&y,conststd::vector<double>&u){intN=static_cast<int>(y.size());std::vector<double>yPred(N,0.0);intstartIdx=std::max(model.na,model.nk+model.nb-1);for(intt=startIdx;t<N;++t){doublepred=0.0;for(inti=1;i<=model.na;++i)pred-=model.a[i-1]*y[t-i];for(inti=0;i<model.nb;++i)pred+=model.b[i]*u[t-model.nk-i];yPred[t]=pred;}returnyPred;}// ------------------------------------------------------------// 示例主程序// ------------------------------------------------------------intmain(){constintN=500;std::vector<double>u(N,0.0),y(N,0.0);std::mt19937rng(42);std::normal_distribution<double>distU(0.0,1.0);std::normal_distribution<double>distNoise(0.0,0.1);for(intt=0;t<N;++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(intt=2;t<N;++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);}autores=estimateARX(y,u,/*na=*/2,/*nb=*/2,/*nk=*/1);if(res.success){std::cout<<"=== 估计结果 ===\n";std::cout<<"a = [";for(doublev:res.model.a)std::cout<<v<<" ";std::cout<<"]\n";std::cout<<"b = [";for(doublev:res.model.b)std::cout<<v<<" ";std::cout<<"]\n";std::cout<<"拟合度: "<<res.fitPercent<<"%\n";std::cout<<"损失函数: "<<res.lossFunction<<"\n";}else{std::cerr<<"估计失败: "<<res.message<<"\n";}return0;}