稀疏技术——使用MATLAB编写

导读

本文源于武老师电力系统仿真课程的其中一个作业。
了解 SPICE的同学可能知道SPICE（SimulationProgram with Integrated Circuit Emphasis）使用的六项技术[改进的节点分析法（ModifiedNodal Analysis），稀疏矩阵解法（Sparse Matrix Solver），牛顿-拉夫逊迭代（Newton-RaphsonIteration），隐性数值积分（ImplicitNumerical Integration），动态步长的瞬态分析（Dynamic Time Step Control），局部截断误差（LocalTruncation Error）]奠定了当代电路仿真软件的基石，而武老师对这六个方面都有详细介绍，所以个人认为即使志不在电力系统，武老师的课程也是个极为有用的课程，对我们理解所有电路仿真软件都有作用。
程序采用MATLAB进行编写。详细的原理性介绍这里不再赘述，可以参考这位师兄的文章：稀疏技术。

核心思想

在这里稀疏矩阵采用散居存储。笔者的程序中核心思想有两个：

1. 将散居存储得到的三个数组合并于一个矩阵，如矩阵第一行存储非零元素的值，第二行存储非零元素在原矩阵中的行号，第三行存储非零元素在原矩阵中的列号。
2. MATLAB中find函数的应用。如定义一个有三行的矩阵Matrix，那么代码
   [xp,yp]=find(Matrix(2,:)==p);
   可以找到Matrix第二行中所有值等于p的元素，并且返回这些元素在Matrix中的行、列号并且存储在数组xp和yp中。

具体代码

1. 对数据采用散居存储
   考虑到电力系统导纳矩阵的对称性，只存储对角线和上三角元素即可，并且将非零元素的行列号以及值存储在一个三行矩阵Matrix中。
2. 因子表分解
   首先寻找Matrix(2,:)中存储着的第i行的非零元素，并且返回该非零元素在Matirx的位置存储在yr数组中；
   接下俩进行因子表分解：
   首先进行规格化计算：根据存储着非零元素在Matrix矩阵第二行的位置的yr数组，可以找到第i行的所有非零元素，然后将非对角线元素值除以该行的对角线线元素值。
   然后接着进行消去计算：
   先对末端节点自边值进行消去计算（最根本的思想就是利用行号=列号这个条件寻找到末端节点的自边值）。
   
   然后对互边进行消去计算。
   
   如果两个末端节点之间不存在边，则需要新增边，即会新出现非零元素，将非零元素直接插到Matrix最后即可。
3. 前代计算
   
   根据图片中的通式编写程序即可。
4. 规格化计算
   
   同样根据图片中的通式编写程序即可。
5. 回代计算
   
   同样根据图片中的通式编写程序即可。

完整代码如下所示：

%% M为传入的矩阵
%  M = [[2,0,0,0,0,0,1,0,0,0,0,1]
%       [0,2,0,0,1,0,0,0,1,0,0,0]
%       [0,0,2,0,0,1,0,0,1,0,0,0]
%       [0,0,0,2,0,0,1,0,0,1,0,0]
%       [0,1,0,0,2,0,0,0,0,0,0,1]
%       [0,0,1,0,0,2,0,0,0,1,0,0]
%       [1,0,0,1,0,0,2,0,0,0,0,0]
%       [0,0,0,0,0,0,0,2,1,0,1,0]
%       [0,1,1,0,0,0,0,1,2,0,0,0]
%       [0,0,0,1,0,1,0,0,0,2,1,0]
%       [0,0,0,0,0,0,0,1,0,1,2,1]
%       [1,0,0,0,1,0,0,0,0,0,1,2]];
%  B = [4;4;4;4;4;4;4;4;5;5;5;5];
 
function Task1(M,B)
[row,column]=size(M); % 获取矩阵M的行数row，列数column
fprintf('注意结算结果是基于对称矩阵情况下计算的\n')
%% 数据采用散居存储，考虑到矩阵对称性，只存储对角线和上三角元素
row_num=[]; % 用来存放非零元素的行号
column_num=[]; % 用来存放非零元素的列号
value=[]; % 用来存放非零元素的值
for i=1:row
    for j=i:column
        if(M(i,j)~=0) 
           value=[value M(i,j)];
           row_num=[row_num i];
           column_num=[column_num j];
        end
    end
end
M=[]; % 释放矩阵M占用的存储空间
% 三者组合成一个矩阵
% 第1行为非零元素值，第2行为非零元素行号，第3行为非零元素列号
Matrix=[value;row_num;column_num];
value=[];row_num=[];column_num=[]; % 释放内存
%% 因子表分解
for i=1:row
    % 输出行数=i的元素在row_num中的位置，y_r为列号
    [~,yr]=find(Matrix(2,:)==i); 
    c_yr=length(yr); % 得到yr数组的大小
    % 规格化计算
    if c_yr>=2
        for j=2:c_yr
            Matrix(1,yr(j))=Matrix(1,yr(j))/Matrix(1,yr(1));
        end
        % 自边的消去计算
        for j=2:c_yr
           % 寻找末端节点自边值的位置
           % 行号和列号相等的点即为末端节点自边值
           p=Matrix(3,yr(j)); % 得到列号
           [~,yp]=find(Matrix(2,:)==p);
           % 末端节点自边值的消去计算
           len_yp=length(yp);
           for b=1:len_yp
               if (Matrix(3,yp(b))==p)
                   Matrix(1,yp(b))=Matrix(1,yp(b))-Matrix(1,yr(j))^2*Matrix(1,yr(1));
               end
           end
        end
        % 互边的消去计算
        if c_yr>=3 % 末端节点在2个以上时，末端节点之间的互边才需要改变或者新增
            for j=2:c_yr
                for k=1:(c_yr-j)
                    % 依次取末端节点任意两条边的列号
                    p1=Matrix(3,yr(j)); 
                    p2=Matrix(3,yr(j+k)); 
                    % 保证边的方向都是编号小的指向编号大的
                    if (p1>=p2) 
                       temp=p1;
                       p1=p2;
                       p2=temp;
                    end
                    % 先查找M(p1,p2)是不是非零元素
                    % 若是非零元素，则直接在Matrix中修改值的大小
                    % 若不是非零元素，则需要新增非零元素在Matrix矩阵中
                    [~,yp1]=find(Matrix(2,:)==p1); % 找到行号为p1的所有非零元素
                    len_yp1=length(yp1); % 得到数组yp1的长度
                    flag=0; % 用来指示两个末端节点之间是否存在边
                    for a=1:len_yp1
                        if (Matrix(3,yp1(a))==p2) 
                           % 两个末端节点之间存在边，则进行修改
                           Matrix(1,yp1(a))=Matrix(1,yp1(a))-Matrix(1,yr(j))*Matrix(1,yr(j+k))*Matrix(1,yr(1));
                           flag=1; % 两个末端节点之间存在边
                           break
                        end
                    end
                    if (flag==0) % 两个末端节点之间不存在边，会新增边
                        temp2=0-Matrix(1,yr(j))*Matrix(1,yr(j+k))*Matrix(1,yr(1));
                        len=length(Matrix(1,:));
                        Matrix(1,len+1)=temp2;
                        Matrix(2,len+1)=p1;
                        Matrix(3,len+1)=p2;
                    end
                end
            end
        end
    end
end
% 得到因子分解表矩阵
D=eye(row);
U=eye(row);
L=eye(row);
temp1=[];
for i=1:length(Matrix(1,:))
    if (Matrix(2,i)==Matrix(3,i))
        temp1=[temp1 Matrix(1,i)];
    else
        U(Matrix(2,i),Matrix(3,i))=Matrix(1,i);
        L(Matrix(3,i),Matrix(2,i))=Matrix(1,i);
    end 
end
for i=1:row
    D(i,i)=temp1(i); 
end
% 对因子表结果进行显示
fprintf('D=');disp(D);
fprintf('U=');disp(U);
%% 前代计算
z=ones(row,1);
temp3=0;
for i=1:row
    for j=1:(i-1)
        [~,yj]=find(Matrix(2,:)==j); % 找到第j行的所有非零元素
        for k=1:length(yj)
            if (Matrix(3,yj(k))==i) % Matrix中存在坐标为(j,i)的非零元素，则进行计算，否则不计算
                temp3=temp3+Matrix(1,yj(k))*z(j,1); 
            end
        end
    end
    z(i,1)=B(i,1)-temp3;
    temp3=0;
end
%% 规格化计算
y=ones(row,1);
for i=1:row
    [~,yi]=find(Matrix(2,:)==i); % 找到第i行的所有非零元素
    for j=1:length(yi)
        if (Matrix(3,yi(j))==i) % 找到对角元元素
           y(i,1)=z(i,1)/Matrix(1,yi(j)); 
        end
    end
end
%% 回代计算
x=ones(row,1);
temp3=0;
for i=1:row
    j=row-i+1;
    for k=(j+1):row
        [~,yj]=find(Matrix(2,:)==j); % 找到第j行的所有非零元素
        for g=1:length(yj)
            if (Matrix(3,yj(g))==k)
                temp3=temp3+Matrix(1,yj(g))*x(k);
            end
        end
    end
    x(j,1)=y(j,1)-temp3;
    temp3=0;
end
fprintf('解x的结果为x=\n');disp(x);
end

结果