线性代数案例

线性代数案例Cayler-Hamilton定理【实验目的】1理解特征多项式的概念 2掌握Cayler-Hamilton定理 【实验要求】掌握生成Vandermonde矩阵的vander命令、求矩阵特征多项式系数的poly()命令、求矩阵范数的norm命令及矩阵多项式运算的polyvalm命令【实验内容】Cayler-Hamilton定理是矩阵理论中的一个比较重要的定理，其内容为：若矩阵A的特征多项式为则有亦即假设矩阵A为Vandermonde矩阵，试验证其满足Cayler-Hamilton定理。【实验方案】Matlab提供了求取矩阵特征多项式系数的函数poly()，但是poly()函数会产生一定的误差，而该误差在矩阵多项式求解中可能导致了巨大的误差，从而得出错误的结论。在实际应用中还有其他简单的数值方法可以精确地求出矩阵的特征多项式系数。例如，下面给出的Fadeev-Fadeeva递推算法也可以求出矩阵的特征多项式。该算法首先给出一个单位矩阵I，并将之赋给，然后对每个k的值分别求出特征多项式参数，并更新矩阵，最终得出矩阵的特征多项式的系数。该算法可以直接由下面的Matlab语句编写一个函数实现：Function c=poly1(A) nr,nc=size(A); if nc=nr % 给出若为方阵，则用Fadeev-Fadeeva算法求特征多项式 I=eye(nc); R=I; c=1 zeros(1,nc); for k=1:nc,c(k+1)=-1/k*trace(A*R);r=A*R+c(k+1)*I; end elseif (nr=1 nc=1) % 给出为向量时，构造矩阵 A=A(isfinite(A);n=length(A) ; % 出去非数或无界的特征根 c=1 zeros(1,n); for j=1:n c(2:(j+1)=c(2:(j+1)-A(j).*c(1:j);endelse % 参数有误则给出错误信息 error (Argument must be a vector or a square matrix.)end.【实验过程】 A = vander(1 2 3 4 5 6 7);运行结果：A = 1 1 1 1 1 1 1 64 32 16 8 4 2 1 729 243 81 27 9 3 1 4096 1024 256 64 16 4 1 15625 3125 625 125 25 5 1 46656 7776 1296 216 36 6 1 117649 16807 2401 343 49 7 1 A运行结果：aa1 = 1.0e+009 * 0.0000 -0.0000 -0.0002 0.0287 1.1589 -6.2505 -2.4223 0.0249如调用新的poly1()函数，则可以得出如下的精确结果。 aa1=poly1(A);b1=polyvalm(aa1,A);norm(B1)运行结果：ans = 0可见，由此得出的B矩阵就会完全等于0，故该矩阵满足Cayley-Hamilton定理。小行星轨道问题【实验目的】1. 掌握线性方程组求解2. 加深对正交变换的理解3. 掌握Matlab软件中的ezplot、zplot命令的区别和适用范围【实验要求】掌握绘制隐函数曲线ezplot命令和彗星状轨迹图comet命令【实验内容】天文学家要确定一颗小行星绕太阳运行的轨道，在轨道平面内建立以太阳为原点的直角坐标系，在两坐标轴上取天文测量单位（一天文单位为地球到太阳的平均距离：9300万里）。在五个不同的时间点对小行星作了观察，测得轨道上五个点的坐标数据如下: 表 2-1 小行星观测数据 x4.55965.08165.55465.96366.2756y0.81451.36851.98952.69253.5265由开普勒第一定律知，小行星轨道为一椭圆。设方程为试确定椭圆的方程并在轨道的平面内以太阳为原点绘出椭圆曲线。并应用坐标平移变换和正交变换将上例题中的二次曲线方程化为标准方程，绘椭圆轨道图，完成小行星运行的动态模拟。【实验方案】（1）二次曲线方程中有五个待定系数：，。将观察所得的五个点坐标数据，代入二次曲线方程得到关于，的线性方程组 求解该方程组得椭圆方程的系数：， 。（2）将椭圆的一般方程写成矩阵形式通过变量变换（平移变换和旋转变换）化为椭圆标准方程。首先化去一次项，然后将二次型化为标准型。为了用平移变换消去一次项，令，（，待定），代入方程整理，得其中，。要化简消去一次项，只须选择，使满足二阶线性方程组将，代入椭圆的一般方程，得令求出特征值极其对应的特征向量。可以取与等价的正交单位向量。构造正交矩阵，利用正交变换得椭圆的标准方程：。椭圆长半轴和短半轴分别为，。【实验过程】(1) MATLAB程序如下：x=4.5596;5.0816;5.5546;5.9636;6.2756;y=0.8145;1.3685;1.9895;2.6925;3.5265;A=x.2,2*x.*y,y.2,2*x,2*y;b=-1;1;1;1;1;a=Ab;syms x y a1 a2 a3 a4 a5fun=a1*x2+2*a2*x*y+a3*y2+2*a4*x+2*a5*y+1;fun=subs(fun,a1,a(1);fun=subs(fun,a2,a(2);fun=subs(fun,a3,a(3);fun=subs(fun,a4,a(4);fun=subs(fun,a5,a(5);ezplot(fun,-1.4,7,-1.5,6.5)运行结果：a=-0.33780.1892-0.38180.46090.4104结果表明：二次曲线方程中的各项系数为=-0.3378，=0.1892，=-0.3818，=0.4609，=0.4104。 图2-2小行星绕太阳运行的轨道(2) MATLAB程序如下：x=4.5596;5.0816;5.5546;5.9636;6.2756;y=0.8145;1.3685;1.9895;2.6925;3.5265;A=x.2,2*x.*y,y.2,2*x,2*y;b=-1;1;1;1;1;ak=Ab;C=ak(1),ak(2);ak(2),ak(3);X=-Cak(4);ak(5);x0=X(1);y0=X(2);X=X;1;D=ak(1),ak(2),ak(4);ak(2),ak(3),ak(5);ak(4),ak(5),1;F=X*D*X;U d=eig(C);a=sqrt(-F/d(1,1);b=sqrt(-F/d(2,2);t=2*pi*(0:5000)/5000;u=a*cos(t);v=b*sin(t);V=U*u;v;x1=V(1,:)+x0;y1=V(2,:)+y0;plot(x1,y1,x,y,*,x0,y0,rO),hold onx2=x1,x1,x1;y2=y1,y1,y1;comet(x2,y2)disp(x0,y0)disp(a,b)图2-3 椭圆轨道图运行结果：2.7213 2.42342.4299 4.3799。结果表明：椭圆标准方程为：。矩阵相似变换在控制理论中的应用【实验目的】1.掌握矩阵的相似变换2.利用矩阵相似变换方法，将控制理论中一般的状态方程变换成某种特殊的形式，以便于更好地进行系统的性质分析3.掌握控制系统的可控标准型、可观察标准型和Jordan标准型【实验要求】掌握Matlab软件中有关相似变换的命令【实验内容】给出系统的相似变换的概念，介绍基于矩阵相似变换的各种标准及变换方法，并用MATLAB编程实现。试求出下面系统的可控标准型:+，并求该状态方程模型的可观测标准型以及Jordan标准型。【实验方案】1. 线性系统的相似变换 假设存在一个非奇异矩阵，且定义了 一个新的状态变量使得，这样关于新状态变量的状态方程模型可以写成 ，且 式中，。在矩阵下的状态变换称为相似变换，称为相似变换矩阵。 2. 单变量控制系统的可控、可观测标准型转换对单变量系统（1）来说，若系统的特征多项式可以写成+ 则可以够造出变换矩阵 这样就可以将原来系统变换成可控制标准型。可以用容易地写出变换矩阵 ；求特征多项式系数，建议用取代 3. 控制系统的Jordan标准型转换 系统的Jordan标准型可以由函数直接求出。值得指出的是，若系统的矩阵含有复数特征值，则用函数不能得出正确结果，应该结合前面Jordan变换的方法手工构造变换矩阵，得出合适的变换系统。【实验过程】（1）得出可控标准型的MATLAB程序 A = -4, -3, 0, -1; -3 , -7, -1, -3; 0, -1, -13, -1; -1, -3, -1, -10; B = 0; -14; 7; 16; C = 0, 0, 12, 0; v = poly1(A);Tc = fliplr(ctrb(A, B) * flipud(hankel(v(end-1: -1: 1)Gc = ss(inv(Tc) * A * Tc, inv(Tc) * B, C * Tc, 0) 变换矩阵和标准型分别为，可观测系统标准型是可控标准型的对偶形式。可观测标准型的变换矩阵为类似于前面的可控标准型变换矩阵，可以由下面语句定义出变换矩阵 % 求特征多项式系数 A=-4, -3, 0, -1; -3, -7, -1, -3; 0, -1, -13, -1; -1, -3, -1, -10 运行结果：A = -4 -3 0 -1 -3 -7 -1 -3 0 -1 -13 -1 -1 -3 -1 -10 B=0, -14, 7, 16 运行结果：B = 0 -14 7 16 C=0, 0, 12, 0 运行结果：C = 0 0 12 0 v = poly1(A) 运行结果：V = 1 34 390 1693 2105 Tc = fliplr(ctrb(A, B) * flipud(hankel(v(end-1: -1: 1) 运行结果：Tc = 5445 776 26 0 -9060 -3909 -433 -14 1415 859 145 7 5400 3178 419 16 Gc = ss(inv(Tc) * A * Tc, inv(Tc) * B, C * Tc, 0) 运行结果：a = x1 x2 x3 x4 x1 0 1 3.553e-015 1.11e-016 x2 0 -4.547e-013 1 0 x3 0 0 4.547e-013 1 x4 -2105 -1693 -390 -34 b = u1 x1 0 x2 0 x3 0 x4 1 c = x1 x2 x3 x4 y1 1.698e+004 1.031e+004 1740 84 d = u1 y1 0 Continuous-time model.（2） 直接求可观测标准型的MATLAB程序 A = -4, -3, 0, -1; -3 , -7, -1, -3; 0, -1, -13, -1; -1, -3, -1, -10; B = 0; -14; 7; 16; C = 0, 0, 12, 0; v = poly1(A);To = inv(fliplr(hankel(v(end-1: -1: 1) * flipud(obsv(A, C)Go = ss(inv(To) * A * To, inv(To) * B, C * To, 0) 变换矩阵和可观测标准型分别为，(3) 求Jordan标准型的MATLAB程序 ； 得出的变换矩阵为 在此变换矩阵下的Jordan标准型为 矩阵的三角分解【实验目的】1.理解矩阵的三角分解（又称为LU分解）2.掌握函数的两种调用方法【实验要求】掌握Matlab软件中有关矩阵LU分解的命令【实验内容】分别用两种方法调用MATLAB中的函数，实现矩阵LU分解问题。【实验方案】矩阵的三角分解又称为LU分解，它的目的是将一个矩阵分解成一个下三角矩阵L和一个上三角矩阵U的乘积，亦即A=LU,其中L和U矩阵可以分别写成，由这两个矩阵可以简单的写出一个矩阵,其中这样产生的矩阵与原来的A矩阵的关系可以写成因此，可以立即得出求取和的递推计算公式及该公式的递推初值为注意，在上述的算法中并未对主元素进行任何选取，因此该算法并不一定数值稳定，因为在运算过程中0或很小的数值可能被用作除数。在MATLAB中也给出了基于主元素的矩阵LU分解函数，该函数的调用格式为分解， 为置换矩阵， 其中，分别为变换后的下三角和上三角矩阵。在MATLAB的函数中考虑了主元素选取的问题，所以该函数一般会给出可靠的结果。又该函数得出的下三角矩阵L并不一定是一个真正的下三角矩阵，因为选取它可能进行了一些元素行的交换，这样主对角线的元素可能不是1，而在矩阵L内存在一个唯一的置换，其各个元素的值均是1.如果想获得有关换行信息，则可以由后一种格式调用函数，这时P为单位阵变换出的置换矩阵，A矩阵可以分解成。【实验过程】 （1）求出三角分解矩阵。 可见，这样得出的矩阵并非下三角矩阵，这是因为再分解过程中采用了主元素交换的方法。现在考虑函数的另一中调用方法。注意，这里得出的P矩阵不是一个单位矩阵，而是单位矩阵的置换矩阵。结合得出的矩阵可以看出，P矩阵的，表明需要将矩阵的第4行换到第2行，表明需要将的第2行换至第3行，将原来第3行换至第4行，这样就可以得出一个真正的下三角矩阵L了。将L,P,U代入并检验，可以精确地还原A矩阵。线性方程组的数值解法【实验目的】1掌握求解上三角形方程组的回代法，求解线性方程组的Gauss消去法、LU分解、Jacobi迭代法、Gauss-Seidel迭代法2比较各种算法的收敛条件及运算效率3使用Matlab编写线性方程组数值解法的M文件及程序【实验要求】1求解上三角形方程组的回代法、求解线性方程组的Gauss消去法、LU分解、Jacobi迭代法、Gauss-Seidel迭代法的原理及步骤2能Matlab编写线性方程组数值解法的M文件及程序【实验内容】用回代法解上三角形方程组，用列主元Gauss消去法、矩阵直接三角分解法、Jacobi迭代法和Gauss-Seidel迭代法解线性方程组。1.求解上三角线性方程组2.用列主元Gauss消去法求解线性方程组3.用矩阵的分解求解线性方程组4.用Jacobi迭代法求解线性方程组5.用Gauss-Seidel迭代法求解线性方程组【实验方案】1. 对于上三角形方程组，其矩阵形式表示为：形如：其中，称为上三角矩阵。若，即，高亦即矩阵是非奇异的，则原方程组有唯一解，且可从最后一个方程解出，即： ， 代入倒数第二个方程，得到： 一般地，设已求得，则由方程的第个方程，可得： 上述求解过种，称为回代过程。回代过程所用乘除法运算次数为，加减法运算次数为。2. Gauss消去法的操作过程如下： 将原线性方程组写成增广矩阵的形式，其中， ， 第一步消元：设，将中的第一行乘以，加到第行上去，可得到同解方程组的增广矩阵：其中， ，这里，符号表示对的第行进行变换，将第一行的倍加到第行上。 第步消元：设，将中的第行对角线以下的元素消为零，即将中的第行乘以，加到第行上去，可得到同解方程组的增广矩阵：其中， ，通常， 称为消元因子，表为主元。 上述做法，直至第步完成，得到同解的方程组：即为上三角矩阵，可表示为式(11.2.2)，这时用回代过程便可求得解向量。而列主元Gauss消去法的基本思想是：设已作步消元，在进行第步消元之前，选出第列中位于对角线及其以下元素绝对值中的最大者，即确定，使得：将的第行和第行互相交换，则元素为新的主元素，其余元素也均以交换后的位置表示。然后再按Gauss消去法进行第步消元。这种方法称为列主元Gauss消去法。此时，消元因子满足：一般能保证舍入误差不扩散，这个方法基本是稳定的。列主元Gauss消去法的运算量除选主元及行交换外和Gauss消去法是相同的。3. Gauss消去过程，实际就是对方程组的增广矩阵连续左乘以得到上三角方程组，即，其中：令，用左乘(11.3.1)式，得，于是，有： 称为单位下三角矩阵，为上三角矩阵。这种矩阵的分解为简单矩阵的乘积，称为矩阵的分解。4. 对于线性方程组设。从第个方程组解出，得到如下同解方程组：建立相应迭代格式：称(上式为Jacobi迭代格式。5. Jacobi迭代格式在逐个求的分量时，当计算到时，分量都已求出，但却被束之高阁，而仍用旧分量进行计算。直观上看，最新算出的分量可能比旧的分量要准确些。因此，设想一旦当新分量已求出，马上就用它来替代，也就是在Jacobi迭代法中求时用代替，这就是Gauss-Seidel迭代法。【实验过程】建立回代法的M函数文件如下：function X=backsub(A,b)% A是一个n阶上三角非奇异阵。% b是一个n维向量。% X是线性方程组AX=b的解。n=length(b);X=zeros(n,1);X(n)=b(n)/A(n,n);for k=n-1:-1:1X(k)=(b(k)-A(k,k+1:n)*X(k+1:n)/A(k,k);end 然后，用上述程序求解上三角线性方程组： 建立一个主程序prog1121.mclcclearA=3,-2,1,-4;0,4,-1,2;0,0,2,3;0,0,0,5;b=8;-3;11;15;backsub(A,b)然后在MATLAB命令窗口运行上述主程序，即： prog1121运行结果：ans = 5 -2 1 3建立列主元 Gauss消去法的函数文件如下：function X=uptrbk(A,b)% A是一个n阶矩阵。% b是一个n维向量。% X是线性方程组AX=b的解。N N=size(A);X=zeros(1,N+1);Aug=A b;for p=1:N-1 Y,j=max(abs(Aug(p:N,p); C=Aug(p,:); Aug(p,:)=Aug(j+p-1,:); Aug(j+p-1,:)=C; if Aug(p,p)=0 A是奇异阵，方程无惟一解 break endfor k=p+1:N m=Aug(k,p)/Aug(p,p); Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1); endend% 这里用前面程序中定义的函数backsub来进行回代。X=backsub(Aug(1:N,1:N)，Aug(1:N,N+1)；再用上述程序求解线性方程组： 建立一个主程序prog1122.mclcclearA=2,4,-6;1,5,3;1,3,2;b=-4;10;5;uptrbk(A,b)然后在MATLAB命令窗口运行上述主程序，即： prog1122运行结果：ans = -3 2 1建立矩阵分解解方程组的函数文件如下：function X=lufact(A,b)% A为n阶矩阵。% b是n维向量。% X是所求的AX=b的解。N,N=size(A);X=zeros(N,1);Y=zeros(N,1);C=zeros(1,N);R=1:N;for p=1:N-1max1,j=max(abs(A(p:N,p);C=A(p,:);A(p,:)=A(j+p-1,:);A(j+p-1,:)=C;d=R(p);R(p)=R(j+p-1);R(j+p-1)=d;if A(p,p)=0 A是奇异阵，方程组无惟一解 breakendfor k=p+1:N mult=A(k,p)/A(p,p); A(k,p)=mult; A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N);endendY(1)=b(R(1);for k=2:N Y(k)=b(R(k)-A(k,1:k-1)*Y(1:k-1);endX(N)=Y(N)/A(N,N);for k=N-1:-1:1 X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N)/A(k,k);end用上述程序求解线性方程组： 建立一个主程序prog1131.mclcclearA=1,3,5,7;2,-1,3,5;0,0,2,5;-2,-6,-3,1;b=1,2,3,4;lufact(A,b)然后在MATLAB命令窗口运行上述主程序，即： prog1131运行结果：ans = 1.3429 0.6857.0000 1.8000 -3还有另一种简便的解法是直接使用MATLAB中对矩阵进行LU分解的命令，然后再用回代法求解原方程组，过程如下：clcclearA=1,3,5,7;2,-1,3,5;0,0,2,5;-2,-6,-3,1;b=1,2,3,4;%对系数矩阵A进行LU分解L,U=lu(A);%用前面程序中定义的函数backsub来进行回代backsub(U,b)建立Jacobi迭代格式的函数文件如下：function X=jacobi(A,b,P,delta,max1)% A是n维非奇异阵。% b是n维向量。% P是初值。% delta是误差界。% max1是给定的迭代最高次数。% X为所求的方程组AX=b的近似解。N=length(b);for k=1:max1for j=1:N X(j)=(b(j)-A(j,1:j-1,j+1:N)*P(1:j-1,j+1:N)/A(j,j); enderr=abs(norm(X-P);P=X； if(err prog1141运行结果：k = 9ans = 3.0000 2.0000 1.0000建立Gauss-Seidel迭代法的M函数文件如下：function X=gseid(A,b,P,delta,max1)% A是n维非奇异阵。% b是n维向量。% P是初值。% delta是误差界。% max1是给定的迭代最高次数。% X为所求的方程组AX=b的近似解。N=length(b);for k=1:max1for j=1:N if j=1 X(1)=(b(1)-A(1,2:N)*P(2:N)/A(1,1); elseif j=N X(N)=(b(N)-A(N,1:N-1)*(X(1:N-1)/A(N,N); else X(j)=(b(j)-A(j,1:j-1)*X(1:j-1)-A(j,j+1:N)*P(j+1:N)/A(j,j); endenderr=abs(norm(X-P);relerr=err/(norm(X)+eps);P=X;if(errdelta)|(relerr prog1142运行结果：err = 8.3073e-005k = 6ans = 2.0000 4.0000 3.0000经过6次迭代得到满足精度要求的解。房屋装修的工资问题【实验目的】1理解矩阵特征值概念 2能根据实际问题，建立模型然后使用Matlab相关命令求解【实验要求】掌握求解特征值的eig命令、生成对角矩阵的diag命令等【实验内容】有三个技术个人分别是木工、电工和管道工，他们准备合作装修自己的新房子。在装修之前约定：每人总共工作20天（包括在自己家）；每人每日的工资平均为100元；每人的日工资应使得每人的总收入和总支出等。需要计算每人的日工资分别是多少，以确定他们的工作日交换是否平衡，如果不平衡，将由谁买单。一个初步的工作日分配方案如下表3-1 工作日分配方案工作日 工种木工电工管道工木工家4212电工家8102管道工家886【实验方案】设木工、电工和管道工的日工资分别为：，。由总收入和总支出相等的约定，建立线性议程组 整理，得 显然问题与矩阵特征值问题有联系，由于矩阵是正矩阵且每列元素之和均为20，所以20是该矩阵的牲值，于是就是属于特征值的特征向量。按约定总工作量决定总工资应该为6000元，则应该有 【实验过程】MATLAB程序如下A=4,2,12;8,10,2;8,8,6; P,D=eig(A); disp(diag(D) II=input(input Index about eigvalu=20:=); if II=0,error(problem have no solution),end alpha=P(:,II); R=alpha./sum(alpha);format bankdaily=300*R pay=A*diag(daily) 运行结果：在MATLAB命令窗口中运行程序，屏幕将显示出A的三个特征值 20.00 -2.00 2.00由于第一个特征值恰好为20，在提示符“input Index about eigvalu=20:=”后输入索引值1。程序继续运行，得出最后计算结果为 daily = 93.94 96.97 109.09 pay = 375.76 193.94 1309.09 751.52 969.70 218.18 751.52 775.76 654.55每人的日工资由变量daily的数据给出。结果表明： 表3-2 日工资列表工种木工电工管道工日工资93.9496.97109.09最后的二维数组给出了二维数组，表明付款明细账，行表示支付，列表示收取。显然第一行相加等于第一列相加，第二行相加等于第二列相加，第三行相加等于第三列相加。 表3-3 工资支付收取方案支付收取木工电工管道工木工375.76193.941309.09电工751.52969.70218.18管道工751.52775.76654.55层次分析法【实验目的】1、掌握正互反矩阵、一致性矩阵等一些基本知识2、理解层次分析法、理想解逼近排序法的基本原理和步骤【实验要求】1、能够对一些定性问题或半定性问题能用层次分析法进行建模求解,并用MATLAB编程实现2、掌握用层次分析法进行建模求解过程Matlab中一些相关命令的用法【实验内容】中国石化股份有限公司西南分公司日前开采的气田主要有孝泉气田、新场气田、合兴场气田等，这些气田最早于1988年投入开采，现已形成较为完善的集输管网系统，日输气能力已达到近400万方/天。川西气田天然气气质好，甲烷含量高达90以上，不含H2S，因此从气井产出的天然气一般只经过采集气站分离、计量、调压后直接外输进入集输管网供给用户。但随着气田开发的进行，气井压力降低，产水逐渐增加，越来越多的气井采用泡沫排水采气工艺进行排水采气，由于站内流程分离器分离效果不佳，加之气井采用泡沫排水采气，部分水分形成泡沫混合在天然气中，导致分离器的分离效果降低，导致天然气含水量增大，对整个集输管网系统运行带来了一系列不利影响。因此，开展川西气田集输管网系统脱水技术应用研究，对提高用户供气质量，降低管网运行压损，消除管内积液对管道的腐蚀具有极为重要的作用。根据不同脱水装置设置点、选址方案以及不同的脱水方法（三甘醇脱水法和硅胶脱水法），得到了6种川西气田集输管网整体脱水方案，见表3-4。表3-4 各个方案的特征指标表指标方案总投资(万元)运营成本(万元/年)处理量（万m3/日）压损（MPa）露点（）保护管线长度（公里）方案13731.49261.904350.2-50方案23726.74256.414350.4-100方案34058.28222.96326.40.2-5398方案44039.75218.64326.40.4-10398方案52785.45222.434000.2-5379方案62800.70183.984000.4-10379由于6个方案都有各自的优点和缺点，如果单纯从技术或经济某一个方面来决定选择一个方案作为最终方案显得不太合理，因此，需要对6个方案进行技术经济综合评价，选择技术经济综合评价相对较高的方案作为最终方案。多方案的技术经济综合评价的主要思路：首先选取综合评价特征指标，由于各特征指标重要性并不完全相同，因而通过层次分析法把各指标重要性从定性转化成定量，作为权重；再通过逼近理想解排序法（TOPSIS法）将各个方案排出优劣次序，确定推荐方案。【实验方案】1 矩阵，矩阵的特征根和特征值；2正互反矩阵、一致性矩阵的定义和性质定义1 若矩阵满足（1），（2）（）则称之为正互反矩阵(易见，)。定义2 若正反矩阵满足 则称A为一致性矩阵。定理1 正互反矩阵的最大特征根必为正实数，其对应特征向量的所有分量均为正实数。的其余特征值的模均严格小于。定理2 若为一致矩阵，则（1）必为正互反矩阵。（2）的转置矩阵也是一致矩阵。（3）的任意两行成比例，比例因子大于零，从而（同样，的任意两列也成比例）。（4）的最大特征值，其中为矩阵的阶。的其余特征根均为零。（5）若的最大特征值对应的特征向量为，则，即定理3 阶正互反矩阵为一致矩阵当且仅当其最大特征根，且当正互反矩阵非一致时，必有。3 层次分析法的基本原理和步骤：层次分析法的基本思路是：将研究对象分解为不同的组成因素,按各因素之间的隶属关系,把它们排成从高到低的若干层次,建立递阶层次结构。对同层的各元素进行两两比较,就每一层次的相对重要性予以定量表示,并利用数学方法确定出每一层次各项因素的权值。层次分析法的流程图如图3-5所示，步骤的详细描述请参见文献数学模型（第三版）姜启源，谢金星。建立层次结构模型构造判断矩阵计算特征值和特征向量一致性检验通过吗？确定指标权重NY图3-5 三标度层次分析法流程图4 理想解逼近排序法（TOPSIS法）的基本原理和步骤：理想解逼近排序法的基本原理为：在维空间中，将方案集中的各个备选方案与理想方案和负理想方案的距离进行比较，既靠近理想方案又远离负理想方案的方案是方案集中的最佳方案；并可以据此排定方案集中各备选方案的优先序。TOPSIS法的流程图如图3-6所示，步骤的详细描述请参阅相关文献。计算规范决策矩阵构成加权规范阵X=xij计算各个方案到理想方案和到负理想方案的距离计算各方案的排序指标值得到各方案的排序图3-6 理想解逼近排序法流程图【实验过程】1.建立层次分析法模型(1) 建立层次结构模型：根据专家意见，选定方案总投资、年运营成本、方案总脱水处理量、脱水前后压损、脱水后干气露点值、保护管线总长为综合评价指标，采用三标度层次分析法构成层次分析模型，对方案进行技术经济评价。综合评价层次结构图3-7：综合评价总目标经济水平技术水平总投资年运营成本处理量压损露点保护管线总长方案1方案6方案2方案3方案4方案5图3-7 层次结构图(2)构造判断矩阵：从层次结构模型的第2层开始，对于从属于（或影响及）上一层每个因素的同一层诸因素，用成对比较法构造判断矩阵。在进行成对比较时，可采用Saaty等人提出的1-9比较尺度。但这里采用了一种较为简化的构造方法：确定各指标综合按重要性大小排列顺序如下：总投资 处理量运营成本 压损露点保护管线长度即，按照层次结构图特征指标顺序，设为第个特征指标（），构造出判断矩阵为（其中特征指标比大则为2，相等则为1，小则为0.5），即：(3)求特征值和特征向量：直接通过计算机编程求解矩阵的特征值和特征向量，找出最大的特征值和其对应的特征向量，最大特征值为 ，其对应的特征向量为： (0.6283，0.3958，0.4986，0.3141，0.2493，0.1979)T计算矩阵特征值和特征向量的matlab命令： A=1,2,2,2,2,2;0.5,1,0.5,2,2,2;0.5,2,1,2,2,2;0.5,0.5,0.5,1,2,2;0.5,0.5,0.5,0.5,1,2;0.5,0.5,0.5,0.5,0.5,1; V,D=eig(A);B=V(:,1)(4)一致性检验：一致性指标的计算公式为：，其中为最大特征根，为矩阵的阶。计算得出，故。根据由Saaty等人给出的层次分析一致性表查得，维数为6的随机一致性指标为1.24，所以 ，检验通过，认为判断矩阵一致性可以接受。其权重为对应向量，将权重化为标准权重，即为：0.2751，0.1733，0.2183，0.1375，0.1092，0.0866 可得到标准权重如表3-8:表3-8 特征指标的标准权重表总投资运营成本处理量压损露点降保护管线长度0.27510.17330.21830.13750.10920.0866计算标准权重的matlab命令： B=B/sum(B)2 理想解逼近排序法（TOPSIS法）（1）计算规范决策矩阵：首先对原始数据（表3-1）标准化处理，得到无量纲数据；然后用向量规范化的方法求得规范决策矩阵。标准化处理采用： 经过以上公式标准化处理后得到矩阵在matlab中编写标准化处理的m文件：clccleary=3731.49 261.90 435 0.2 -5 0; 3726.74 256.41 435 0.4 -10 0; 4058.28 222.96 326.4 0.2 -5 398; 4039.75 218.64 326.4 0.4 -10 398; 2785.45 222.43 400 0.2 -5 379; 2800.70 183.98 400 0.4 -10 379;c=zeros(6,6);for j=1:6 for i=1:6 c(i,j)=y(i,j)/sqrt(sum(y(:,j).2); endendc（2）构成加权规范阵。设由层次分析法得到的权重为则，其中， 即得到为：计算的matlab程序命令：for i=1:6