C++代码实现MATLAB中的rss函数功能

// rss.cpp ------ MATLAB rss 的纯 C 实现无第三方库// 编译: g -stdc17 -O2 rss.cpp -o rss#includeiostream#includevector#includerandom#includecmath#includeiomanip#includealgorithm#includestringusingMatrixstd::vectorstd::vectordouble;// -------- 基础矩阵运算 --------Matrixrandom_matrix(introws,intcols,std::mt19937gen){std::normal_distributiondoubledist(0.0,1.0);MatrixM(rows,std::vectordouble(cols));for(inti0;irows;i)for(intj0;jcols;j)M[i][j]dist(gen);returnM;}Matrixmatmul(constMatrixA,constMatrixB){intnA.size(),mB[0].size(),kB.size();MatrixC(n,std::vectordouble(m,0.0));for(inti0;in;i)for(intj0;jm;j)for(intl0;lk;l)C[i][j]A[i][l]*B[l][j];returnC;}Matrixtranspose(constMatrixA){intnA.size(),mA[0].size();MatrixB(m,std::vectordouble(n));for(inti0;in;i)for(intj0;jm;j)B[j][i]A[i][j];returnB;}Matrixidentity(intn){MatrixI(n,std::vectordouble(n,0.0));for(inti0;in;i)I[i][i]1.0;returnI;}// -------- Householder QR 分解: A Q * R --------voidqr_decompose(constMatrixA,MatrixQ,MatrixR){intnA.size(),mA[0].size();RA;Qidentity(n);for(intk0;kstd::min(n,m);k){doublenorm0;for(intik;in;i)normR[i][k]*R[i][k];normstd::sqrt(norm);if(norm1e-14)continue;doublealpha(R[k][k]0)?-norm:norm;std::vectordoublev(n,0.0);for(intik;in;i)v[i]R[i][k];v[k]-alpha;doublevnorm20;for(intik;in;i)vnorm2v[i]*v[i];if(vnorm21e-14)continue;// 左乘 Householder更新 Rfor(intj0;jm;j){doubledot0;for(intik;in;i)dotv[i]*R[i][j];dot*2.0/vnorm2;for(intik;in;i)R[i][j]-dot*v[i];}// 右乘 Householder累积 Q Q * H_kfor(inti0;in;i){doubledot0;for(intjk;jn;j)dotQ[i][j]*v[j];dot*2.0/vnorm2;for(intjk;jn;j)Q[i][j]-dot*v[j];}}}// 通过随机矩阵 QR 分解生成均匀分布的随机正交矩阵Matrixrandom_orthogonal(intn,std::mt19937gen){Matrix Mrandom_matrix(n,n,gen);Matrix Q,R;qr_decompose(M,Q,R);// 使 R 的对角元为正保证 Q 唯一for(inti0;in;i)if(R[i][i]0)for(intj0;jn;j)Q[j][i]-Q[j][i];returnQ;}// -------- 状态空间结构 --------structStateSpace{Matrix A,B,C,D;};// 等价于 MATLAB: sys rss(n, p, m)StateSpacerss(intn,intp1,intm1,unsignedseed42){std::mt19937gen(seed);std::normal_distributiondoublend(0.0,1.0);std::uniform_real_distributiondoubleud(0.0,1.0);// 1) 构造随机实 Schur 形式 TMatrixT(n,std::vectordouble(n,0.0));std::vectorintblock_type(n,0);// 11x1块; 2,32x2块的两个位置inti0;while(in){if(in-1ud(gen)0.5){// 复共轭极点对 a ± bj (a 0)doublea-std::abs(nd(gen))-0.1;doublebnd(gen);if(std::abs(b)0.1)b0.5;// 避免退化为实极点T[i][i]a;T[i][i1]b;T[i1][i]-b;T[i1][i1]a;block_type[i]2;block_type[i1]3;i2;}else{// 实极点 (负)T[i][i]-std::abs(nd(gen))-0.1;block_type[i]1;i1;}}// 填充严格上三角部分2x2 块内的元素已赋值跳过for(intr0;rn;r)for(intcr1;cn;c){if(block_type[r]2block_type[c]3cr1)continue;T[r][c]nd(gen);}// 2) 随机正交矩阵 QMatrix Qrandom_orthogonal(n,gen);// 3) A Q * T * Q^T ------ 特征值 T 的对角块特征值全部稳定Matrix Amatmul(matmul(Q,T),transpose(Q));// 4) 随机 B, C, D 并归一化Matrix Brandom_matrix(n,m,gen);Matrix Crandom_matrix(p,n,gen);Matrix Drandom_matrix(p,m,gen);for(intj0;jm;j){doublenrm0;for(intk0;kn;k)nrmB[k][j]*B[k][j];nrmstd::sqrt(nrm);if(nrm1e-12)for(intk0;kn;k)B[k][j]/nrm;}for(intj0;jp;j){doublenrm0;for(intk0;kn;k)nrmC[j][k]*C[j][k];nrmstd::sqrt(nrm);if(nrm1e-12)for(intk0;kn;k)C[j][k]/nrm;}return{A,B,C,D};}// -------- 打印 --------voidprint_matrix(constMatrixM,conststd::stringname){std::coutname \n;for(constautorow:M){for(doublev:row)std::coutstd::setw(10)std::setprecision(4)std::fixedv ;std::cout\n;}std::cout\n;}intmain(){std::cout C rss(3, 1, 1) \n\n;StateSpace sysrss(3,1,1,12345);print_matrix(sys.A,A);print_matrix(sys.B,B);print_matrix(sys.C,C);print_matrix(sys.D,D);std::coutA 由 A Q*T*Q^T 构造T 为随机实 Schur 形式\n其对角块1x1 或 2x2的特征值实部均严格为负\n因此 A 保证稳定Hurwitz。\n;return0;}