最短路算法,matlab程序
利用Matlab编程计算最短途径及中位点选址
§19. 利用Matlab 编程计算最短途径及中位点选址一、最短路问题两个指定极点之间的最短途径。
例如，给出了一个连接假设干个城镇的铁路网络，在那个网络的两个指定城镇间，找一条最短铁线路。
以各城镇为图G 的极点，两城镇间的直通铁路为图G 相应两极点间的边，得图G 。
对G 的每一边e ，赋以一个实数)(e w —直通铁路的长度，称为e 的权，取得赋权图G 。
G 的子图的权是指子图的各边的权和。
问题确实是求赋权图G 中指定的两个极点00,v u 间的具最小权的轨。
这条轨叫做00,v u 间的最短路，它的权叫做00,v u 间的距离，亦记作),(00v u d 。
求最短路已有成熟的算法：迪克斯特拉（Dijkstra ）算法，其大体思想是按距0u 从近到远为顺序，依次求得0u 到G 的各极点的最短路和距离，直至0v （或直至G 的所有极点），算法终止。
为幸免重复并保留每一步的计算信息，采纳了标号算法。
下面是该算法。
(i) 令0)(0=u l ，对0u v ≠，令∞=)(v l ，}{00u S =，0=i 。
(ii) 对每一个i S v ∈（i i S V S \=），用)}()(),({min uv w u l v l iS u +∈代替)(v l 。
计算)}({min v l iS v ∈，把达到那个最小值的一个极点记为1+i u ，令}{11++=i i i u S S 。
(iii). 若1||-=V i ，停止；假设1||-<V i ，用1+i 代替i ，转(ii)。
算法终止时，从0u 到各极点v 的距离由v 的最后一次的标号)(v l 给出。
在v 进入i S 之前的标号)(v l 叫T 标号，v 进入i S 时的标号)(v l 叫P 标号。
算法确实是不断修改各项点的T 标号，直至取得P 标号。
假设在算法运行进程中，将每一极点取得P 标号所由来的边在图上标明，那么算法终止时，0u 至各项点的最短路也在图上标示出来了。
最短路dijkstra算法Matlab程序
function [c0,c,path0,path]=dijkstra(s,t,C,flag)% Use the Dijkstra's algorithm to find the shortest path from% s to t and can also find the shortest path between s and all% the other points.% Reference: Graph Theory with Applications by J. A. Bondy and% U. S. R. Murty.% Input -- s is the starting point and also is the point s.% -- t is the given terminal point and is the point t.% -- C \in R^{n \times n}is the cost matrix, where% C(i,j)>=0 is the cost from point i to point j.% If there is no direct connection between point i and% j, C(i,j)=inf.% -- flag: if flag=1, the function just reports the% shortest path between s and t; if flag~=1, the% function reports the shortest path between s and t,% and the shortest paths between s and other points.% Output -- c0 is the minimal cost from s to t.% -- path0 denotes the shortest path form s to t.% -- c \in R{1\times n} in which the element i is the% minimal cost from s to point i.% -- path \in R^{n \times n} in which the row i denotes% the shortest path from s to point i.% Copyright by MingHua Xu(徐明华), Changhzou University, 27 Jan. 2014. s=floor(s);t=floor(t);n=size(C,1);if s<1 || t < 1 || s > n || t > nerror(' The starting point and the terminal point exceeds the valid range');endif t==sdisp('The starting point and the terminal point are the same points');endlabel=ones(1,n)*inf;label(s)=0;S=[s];Sbar=[1:s-1,s+1:n];c0=0;path=zeros(n,n);path(:,1)=s;c=ones(1,n)*inf;parent=zeros(1,n);i=1; % number of points in point set S.while i<n% for each point in Sbar, replace label(Sbar(j)) by% min(label(Sbar(j)),label(S(k))+C(S(k),Sbar(j)))for j=1:n-ifor k=1:iif label(Sbar(j)) > label(S(k))+C(S(k),Sbar(j))label(Sbar(j))=label(S(k))+C(S(k),Sbar(j));parent(Sbar(j))=S(k);endendend% Find the minmal label(j), j \in Sbar.temp=label(Sbar(1));son=1;for j=2:n-iif label(Sbar(j))< temptemp=label(Sbar(j));son=j;endend% update the point set S and SbarS=[S,Sbar(son)];Sbar=[Sbar(1:son-1),Sbar(son+1:n-i)];i=i+1;% if flag==1, just output the shortest path between s and t.if flag==1 && S(i)==tson=t;temp_path=[son];if son~=swhile parent(son)~=sson=parent(son);temp_path=[temp_path,son];endtemp_path=[temp_path,s];endtemp_path=fliplr(temp_path);m=size(temp_path,2);path0(1:m)=temp_path;c_temp=0;for j=1:m-1c_temp=c_temp+C(temp_path(j),temp_path(j+1));endc0=c_temp;path(t,1:m)=path0;c(t)=c0;returnendend% Form the output resultsfor i=1:nson=i;temp_path=[son];if son~=swhile parent(son)~=sson=parent(son);temp_path=[temp_path,son];endtemp_path=[temp_path,s];endtemp_path=fliplr(temp_path);m=size(temp_path,2);path(i,1:m)=temp_path;c_temp=0;for j=1:m-1c_temp=c_temp+C(temp_path(j),temp_path(j+1));endc(i)=c_temp;c0=c(t);path0=path(t,:);endreturn。
Floyd最短路算法的MATLAB程序
Floyd最短路算法的MATLAB程序Floyd最短路算法的MATLAB程序2006-08-17 20:14%floyd.m%采用floyd算法计算图a中每对顶点最短路 %d是矩离矩阵%r是路由矩阵function [d,r]=floyd(a)n=size(a,1);d=a;for i=1:nfor j=1:nr(i,j)=j;endendrfor k=1:nfor i=1:nfor j=1:nif d(i,k)+d(k,j)<d(i,j)< p="">d(i,j)=d(i,k)+d(k,j); r(i,j)=r(i,k)endendendkdrendvoid Dijkstral(int v0){int i;bool s[MAX_VEX];for(i=0;i<dim;i++)< p="">{d[v0][i]=map[v0][i];s[i]=false;if((i!=0)&&(d[v0][i]<inf))< p="">p[v0][i]=v0;elsep[v0][i]=-1;}s[v0]=true;d[v0][v0]=0;for(i=0;i<dim;i++)< p="">{double min=INF;int u=v0;for(int j=0;j<dim;j++)< p="">if(!s[j]&&d[v0][j]<min)< p="">{u=j;min=d[v0][j];}s[u]=true;for(int w=0;w<dim;w++)< p="">{if((!s[w])&&(d[v0][w]>d[v0][u]+map[u][w])) {d[v0][w]=d[v0][u]+map[u][w];p[v0][w]=u;}}}}Justin Hou介绍寻找最有价值路径（c语言）描述：从上（入口）往下行走，直到最下节点（出口）结束，将所经节点上的数值相加，要求找到一条最有价值路径（既是路径总数值最大）并输出总数值。
D_i_j_k_s_t_r_a最短路算法MATLAB程序_
从起点sb到终点db通用的Dijkstra标号算法程序function [mydistance,mypath]=mydijkstra(a,sb,db);% 输入：a—邻接矩阵，a（i，j)是指i到j之间的距离，可以是有向的% sb—起点的标号, db—终点的标号% 输出：mydistance—最短路的距离, mypath—最短路的路径n=size(a,1); visited(1:n) = 0;distance(1:n) = inf; distance(sb) = 0; %起点到各顶点距离的初始化visited(sb)=1; u=sb; %u为最新的P标号顶点parent(1:n) = 0; %前驱顶点的初始化for i = 1: n-1id=find(visited==0); %查找未标号的顶点for v = idif a(u, v) + distance(u) < distance(v)distance(v) = distance(u) + a(u, v); %修改标号值parent(v) = u;endendtemp=distance;temp(visited==1)=inf; %已标号点的距离换成无穷[t, u] = min(temp); %找标号值最小的顶点visited(u) = 1; %标记已经标号的顶点endmypath = [];if parent(db) ~= 0 %如果存在路!t = db; mypath = [db];while t ~= sb %从终点db开始回溯p = parent(t);mypath = [p mypath];t = p;endendmydistance = distance(db);例题：运筹学教材P205 第七题D=[0 3 6 1 inf inf inf inf;inf 0 2 inf 4 6 inf inf;inf inf 0 inf inf 5 inf inf;inf inf 4 0 inf 3 6 inf;inf inf inf inf 0 inf inf 7;inf inf inf inf inf 0 7 11;inf inf inf inf inf inf 0 8;inf inf inf inf inf inf inf 0]在command window输入：[mydistance,mypath]=mydijkstra(D,1,8);。
最短距离算法matlab
最短距离算法matlab最短距离算法在计算机科学领域中被广泛应用于寻找两个点之间最短路径的问题。
无论是在社交网络分析、物流规划还是路线规划等应用中，最短距离算法都扮演着核心角色。
本文将介绍最短距离算法在MATLAB中的实现，并详细说明算法的原理和步骤。
1. 引言最短距离算法（Shortest Path Algorithm）用于计算两个点之间路径的最短距离。
该算法通过对图模型进行分析和搜索，找到连接两个点的最短路径。
最短路径问题可以通过多种算法求解，例如迪杰斯特拉算法（Dijkstra's Algorithm）和弗洛伊德–沃尓沃兹算法（Floyd-Warshall Algorithm）等。
在本文中，我们将着重解释迪杰斯特拉算法的MATLAB 实现。
2. 算法原理迪杰斯特拉算法是一种用于解决最短路径问题的贪心算法。
它通过逐步构建最短路径树来确定最短路径。
算法的核心思想是通过不断更新节点的距离和前驱节点来找到最短路径。
3. 算法步骤以下是迪杰斯特拉算法的步骤：步骤1: 创建一个空的距离和前驱节点的集合。
距离集合用于存储源节点到其他节点的距离，前驱节点集合用于存储最短路径中每个节点的前驱节点。
距离集合初始化为无穷大，源节点的距离初始化为0。
前驱节点集合初始化为空。
步骤2: 选择源节点，并将其添加到已访问的节点集合中。
步骤3: 对于源节点的邻居节点，更新其距离和前驱节点。
如果新的距离小于当前距离，则更新距离和前驱节点。
步骤4: 从距离集合中选择最小距离的节点，并将其添加到已访问的节点集合中。
该节点将成为当前节点。
步骤5: 重复步骤3和步骤4，直到所有节点都被访问。
步骤6: 根据前驱节点集合构建最短路径。
4. MATLAB代码实现下面是迪杰斯特拉算法在MATLAB中的实现：MATLABfunction [dist, path] = dijkstra(graph, source_node,destination_node)n = size(graph, 1); 图中节点的数量dist = inf(1, n); 节点到源节点的距离dist(source_node) = 0; 源节点到自身的距离为0visited = false(1, n); 记录节点是否已被访问path = zeros(1, n); 记录最短路径的前驱节点for i = 1:n-1[min_dist, current_node] = min(dist); 选择距离最小的节点visited(current_node) = true; 标记该节点为已访问for j = 1:nif ~visited(j) && graph(current_node, j) > 0 如果节点未被访问且节点之间存在连接new_dist = dist(current_node) +graph(current_node, j); 更新距离if new_dist < dist(j)dist(j) = new_dist; 更新距离集合中的节点距离path(j) = current_node; 更新前驱节点endendendend构建最短路径path_len = 1;current_node = destination_node;while current_node ~= source_nodepath_len = path_len + 1;current_node = path(current_node);endpath = zeros(1, path_len);current_node = destination_node;for i = path_len:-1:1path(i) = current_node;current_node = path(current_node);endend5. 实例演示现在我们将使用上述实现的迪杰斯特拉算法来计算一个简单图中两个节点之间的最短路径。
metlab最短路标号法代码及结果
metlab最短路标号法代码及结果MATLAB 最短路-标号法代码函数代码function [min,path]=dijkstra(w,start,terminal)n=size(w,1); label(start)=0; f(start)=start;for i=1:nif i~=startlabel(i)=inf;end, ends(1)=start; u=start;while length(s)for i=1:nins=0;for j=1:length(s)if i==s(j)ins=1;end, endif ins==0v=i;if label(v)>(label(u)+w(u,v))label(v)=(label(u)+w(u,v)); f(v)=u;end, end, endv1=0;k=inf;for i=1:nins=0;for j=1:length(s)if i==s(j)ins=1;end, endif ins==0v=i;if k>label(v)k=label(v); v1=v;u=v1;endmin=label(terminal); path(1)=terminal; i=1;while path(i)~=startpath(i+1)=f(path(i));i=i+1 ;endpath(i)=start;L=length(path);path=path(L:-1:1);脚本代码weight=[0 3 5 inf inf inf inf;3 0 inf 2 inf inf inf ;5 inf 0 1.5 4 inf inf;inf 2 1.5 0 inf 3 inf;inf inf 4 inf 0 0.5 3;inf inf inf 3 0.5 0 1.5;inf inf inf inf 3 1.5 inf;]; [dis, path]=dijkstra(weight, i,j) 输出结果>> weight=[0 3 5 inf inf inf inf;3 0 inf 2 inf inf inf ;5 inf 0 1.5 4 inf inf;inf 2 1.5 0 inf 3 inf;inf inf 4 inf 0 0.5 3;inf inf inf 3 0.5 0 1.5;inf inf inf inf 3 1.5 inf;]; [dis, path]=dijkstra(weight,1,4)dis =5path =0 3 5 inf inf inf inf;3 0 inf 2 inf inf inf ;5 inf 0 1.5 4 inf inf;inf 2 1.5 0 inf 3 inf;inf inf 4 inf 0 0.5 3;inf inf inf 3 0.5 0 1.5;inf inf inf inf 3 1.5 inf;]; [dis, path]=dijkstra(weight,1,5) dis = 8.5000path =1 2 4 6 5>> weight=[0 3 5 inf inf inf inf;3 0 inf 2 inf inf inf ;5 inf 0 1.5 4 inf inf;inf 2 1.5 0 inf 3 inf;inf inf 4 inf 0 0.5 3;inf inf inf 3 0.5 0 1.5;inf inf inf inf 3 1.5 inf;]; [dis, path]=dijkstra(weight,1,7)dis =9.5000path =1 2 4 6 7>> weight=[0 3 5 inf inf inf inf;3 0 inf 2 inf inf inf ;5 inf 0 1.5 4 inf inf;inf 2 1.5 0 inf 3 inf;inf inf 4 inf 0 0.5 3;inf inf inf 3 0.5 0 1.5;inf inf inf inf 3 1.5 inf;]; [dis, path]=dijkstra(weight,4,5) dis =>> weight=[0 3 5 inf inf inf inf;3 0 inf 2 inf inf inf ;5 inf 0 1.5 4 inf inf;inf 2 1.5 0 inf 3 inf;inf inf 4 inf 0 0.5 3;inf inf inf 3 0.5 0 1.5;inf inf inf inf 3 1.5 inf;]; [dis, path]=dijkstra(weight,4,7) dis = 4.5000path =4 6 7>> weight=[0 3 5 inf inf inf inf;3 0 inf 2 inf inf inf ;5 inf 0 1.5 4 inf inf;inf 2 1.5 0 inf 3 inf; inf inf 4 inf 0 0.5 3; inf inf inf 3 0.5 0 1.5;inf inf inf inf 3 1.5 inf;]; [dis, path]=dijkstra(weight,5,7)dis =2path =5 6 7>> weight=[0 3 5 inf inf inf inf;3 0 inf 2 inf inf inf ;5 inf 0 1.5 4 inf inf;inf 2 1.5 0 inf 3 inf;inf inf 4 inf 0 0.5 3;inf inf inf 3 0.5 0 1.5;inf inf inf inf 3 1.5 inf;path =7 6 4 2 1>> weight=[0 3 5 inf inf inf inf;3 0 inf 2 inf inf inf ;5 inf 0 1.5 4 inf inf;inf 2 1.5 0 inf 3 inf;inf inf 4 inf 0 0.5 3;inf inf inf 3 0.5 0 1.5;inf inf inf inf 3 1.5 inf;]; [dis, path]=dijkstra(weight,5,1) dis =8.5000path =5 6 4 2 1>> weight=[0 3 5 inf inf inf inf;3 0 inf 2 inf inf inf ;5 inf 0 1.5 4 inf inf;inf 2 1.5 0 inf 3 inf;inf inf 4 inf 0 0.5 3;inf inf inf 3 0.5 0 1.5;inf inf inf inf 3 1.5 inf;]; [dis, path]=dijkstra(weight,4,1) dis =5path =4 2 1>> weight=[0 3 5 inf inf inf inf;3 0 inf 2 inf inf inf ;5 inf 0 1.5 4 inf inf;inf inf inf 3 0.5 0 1.5;inf inf inf inf 3 1.5 inf;]; [dis, path]=dijkstra(weight,5,4)dis =3.5000path =5 6 4>> weight=[0 3 5 inf inf inf inf;3 0 inf 2 inf inf inf ;5 inf 0 1.5 4 inf inf;inf 2 1.5 0 inf 3 inf;inf inf 4 inf 0 0.5 3;inf inf inf 3 0.5 0 1.5;inf inf inf inf 3 1.5 inf;]; [dis, path]=dijkstra(weight,7,4) dis = 4.5000path =7 6 4>> weight=[0 3 5 inf inf inf inf;3 0 inf 2 inf inf inf ;5 inf 0 1.5 4 inf inf;inf 2 1.5 0 inf 3 inf;inf inf 4 inf 0 0.5 3;inf inf inf 3 0.5 0 1.5;inf inf inf inf 3 1.5 inf;]; [dis, path]=dijkstra(weight,7,5) dis = 2path =7 6 5。
短路计算matlab程序
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Power System Analysis' Project: Matlab Program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear allclose allclca=exp(j*2*pi/3);A=[1 1 1 ;1 a^2 a ;1 a a^2];Sref=900e6;Uref=300e3;Iref=Sref/(sqrt(3)*Uref);% data from the different documents:Positive;Negative;Zero;% run power flow:opt = mpoption('ENFORCE_Q_LIMS', 1);results=runpf(case9gs,opt);Vp=[results.bus(:,8).*(cos(pi/180*results.bus(:,9))+j*sin(pi/180*result s.bus(:,9)))];clc;% Menu:Bucle=1;while (Bucle==1)% Menu:disp('_________________________________________________________________ _______________________')disp('%%%%%% Program to calculate the symetrical and unsymetrical short circuit faults %%%%%%');disp(' 1) Symetrical fault (three phase short circuit fault)');disp(' 2) Unsymetrical fault (A single phase to ground fault)');disp(' 0) For exit the program');obtion=input('Choose the option you want to be calculated by typing its number:','s');disp('_________________________________________________________________ _______________________')% Calculation of the faultswitch (obtion)%%case'1'disp('You Have choose: SYMETRICAL FAULT');disp(' ');disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA*')disp('*******************************************************************************')bus=input ('Please, number the bus where the three phase short circuit fault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number X will have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase. Rf=');% Calculation:If=zeros(8,1);If(bus)=Vp(bus)/(Zpos(bus,bus)+Rf);Vs=Vp-Zpos*If;% Showing the Results:disp(' ');disp('*******************************************************************************')disp('* RESULTS*')disp('*******************************************************************************')FaultCurrent=sprintf('The short circuit currents is: %gA. with angle %g?, Iref*abs(If(bus)),180/pi*angle(If(bus)));disp(FaultCurrent);figure();compass(If(bus),'blue');hold oncompass(If(bus)*a,'red');compass(If(bus)*a^2,'green');text(-abs(If(bus)),-abs(If(bus)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp ('_____________________________________________')disp (' Bus VAngle[Degree]')disp ('_____________________________________________')VoltageTokke= sprintf(' Tokke %gV %g?,Uref*abs(Vs(1)),180/pi*angle(Vs(1)));disp(VoltageTokke) VoltageVinje= sprintf(' Vinje %gV %g?,Uref*abs(Vs(2)),180/pi*angle(Vs(2)));disp(VoltageVinje) VoltageSonga= sprintf(' Songa %gV %g?,Uref*abs(Vs(3)),180/pi*angle(Vs(3)));disp(VoltageSonga) VoltageVemork= sprintf(' Vemork %gV %g?,Uref*abs(Vs(4)),180/pi*angle(Vs(4)));disp(VoltageVemork) VoltageRjukan= sprintf(' Rjunkan %gV %g?,Uref*abs(Vs(5)),180/pi*angle(Vs(5)));disp(V oltageRjukan) VoltageFlesaker=sprintf(' Flesaker %gV %g?,Uref*abs(Vs(6)),180/pi*angle(Vs(6)));disp(VoltageFlesaker) VoltageTveiten= sprintf(' Tveiten %gV %g?,Uref*abs(Vs(7)),180/pi*angle(Vs(7)));disp(Volta geTveiten) VoltageRod= sprintf(' Rod %gV %g?,Uref*abs(Vs(8)),180/pi*angle(Vs(8)));disp(VoltageRod) disp ('_____________________________________________')figure()subplot(2,4,1)compass(Vs(1),'b');hold oncompass(Vs(1)*a,'r');compass(Vs(1)*a^2,'g');text(-abs(Vs(1)),-abs(Vs(1)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(Vs(2));hold oncompass(Vs(2)*a,'r');compass(Vs(2)*a^2,'g');text(-abs(Vs(2)),-abs(Vs(2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(Vs(3));hold oncompass(Vs(3)*a,'r');compass(Vs(3)*a^2,'g');text(-abs(Vs(3)),-abs(Vs(3)),'Songa`s Voltage vectors');subplot(2,4,4)compass(Vs(4));hold oncompass(Vs(4)*a,'r');compass(Vs(4)*a^2,'g');text(-abs(Vs(4)),-abs(Vs(4)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(Vs(5));hold oncompass(Vs(5)*a,'r');compass(Vs(5)*a^2,'g');text(-abs(Vs(5)),-abs(Vs(5)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(Vs(6));hold oncompass(Vs(6)*a,'r');compass(Vs(6)*a^2,'g');text(-abs(Vs(6)),-abs(Vs(6)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(Vs(7));hold oncompass(Vs(7)*a,'r');compass(Vs(7)*a^2,'g');text(-abs(Vs(7)),-abs(Vs(7)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(Vs(8));hold oncompass(Vs(8)*a,'r');compass(Vs(8)*a^2,'g');text(-abs(Vs(8)),-abs(Vs(8)),'Rod`s Voltage vectors');disp(' ')disp('*******************************************************************************')disp('* CURRENTS*')disp('*******************************************************************************')Q=input ('Do you want to know any of the currents through two connected points? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch(k)=(Vs(B1(k))-Vs(B2(k)))*(-Ypos(B1(k),B2(k)));endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('*******************************************************************************')disp('* CURRENTS RESULTS*')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');endfor K=1:kCurrent=sprintf ('From bus no %g to bus no %g the current flowing is: %gA %g?,B1(K),B2(K),Iref*abs(Ibranch(K)),180/pi*angle(Ibranch( K)));disp (Current)figure();compass(Ibranch(K),'blue');hold oncompass(Ibranch(K)*a,'red');compass(Ibranch(K)*a^2,'green');text(-abs(Ibranch(K)),-abs(Ibranch(K)),'Current between selected nodes');end%%case'2'disp('You Have choose: UNSYMETRICAL FAULT');disp('The program only calcutes single line to ground faults, it is consider to happen in phase "a"')disp(' ')disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA*')disp('*******************************************************************************')bus=input ('Please, number the bus where the single phase to ground fault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number X will have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase. Rf=');% Calculation of the currents:Ifphaseground=3*Vp(bus)/(Zpos(bus,bus)+Zneg(bus,bus)+Zzero(bus,bus)+3*R f);%Fase aIfa=zeros(8,1);Ifa(bus)=Ifphaseground;Ifa1=Ifa/3;Ifa2=Ifa/3;Ifa0=Ifa/3;for k=1:8Ifault=A*[Ifa0(k);Ifa1(k);Ifa2(k)];IF(k,:)=Ifault';end% Voltages;%Fase a:V1=Vp-Zpos*Ifa1;V2=-Zneg*Ifa2;V0=-Zzero*Ifa0;for k=1:8U=A*[V0(k);V1(k);V2(k)]V(k,:)=transpose(U)end% Showing the Results:disp(' ');disp('***************************************************************** **************')disp('* RESULTS*')disp('***************************************************************** **************')disp(' ')FaultCurrent=sprintf('The short circuit currents is: %g A. with angle %g ?, Iref*abs(Ifa(bus)),180/pi*angle(Ifa(bus)));disp(FaultCurrent);figure();compass(IF(bus,1),'red');hold oncompass(IF(bus,2),'blue');compass(IF(bus,3),'green');text(-abs(IF(bus,1)),-abs(IF(bus,1)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp('_________________________________________________________________ _________________________________')disp(' Bus Va VbVc');disp('__________________________________________________________________________________________________')VoltageTokke= sprintf('Tokke %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(1,1)),Uref/sqrt(3)*abs(V(1,2) ),180/pi*angle(V(1,2)),Uref/sqrt(3)*abs(V(1,3)),180/pi*angle(V(1,3)));d isp(VoltageTokke)VoltageVinje= sprintf('Vinje %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(2,2) ),180/pi*angle(V(2,2)),Uref/sqrt(3)*abs(V(2,3)),180/pi*angle(V(2,3)));d isp(VoltageVinje)VoltageSonga= sprintf('Songa %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(3,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(3,2) ),180/pi*angle(V(3,2)),Uref/sqrt(3)*abs(V(3,3)),180/pi*angle(V(3,3)));d isp(VoltageSonga)VoltageVemork= sprintf('Vemork %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(4,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(4,2) ),180/pi*angle(V(4,2)),Uref/sqrt(3)*abs(V(4,3)),180/pi*angle(V(4,3)));d isp(VoltageVemork)VoltageRjukan= sprintf('Rjukan %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(5,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(5,2) ),180/pi*angle(V(5,2)),Uref/sqrt(3)*abs(V(5,3)),180/pi*angle(V(5,3)));d isp(VoltageRjukan)VoltageFlesaker=sprintf('Flesaker %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(6,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(6,2) ),180/pi*angle(V(6,2)),Uref/sqrt(3)*abs(V(6,3)),180/pi*angle(V(6,3)));d isp(VoltageFlesaker)VoltageTveiten= sprintf('Tveiten %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(7,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(7,2) ),180/pi*angle(V(7,2)),Uref/sqrt(3)*abs(V(7,3)),180/pi*angle(V(7,3)));d isp(VoltageTveiten)VoltageRod= sprintf('Rod %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(8,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(8,2) ),180/pi*angle(V(8,2)),Uref/sqrt(3)*abs(V(8,3)),180/pi*angle(V(8,3)));d isp(VoltageRod)disp('_________________________________________________________________ _________________________________')figure()subplot(2,4,1)compass(V(1,2),'b');hold oncompass(V(1,1),'r');compass(V(1,3),'g');text(-abs(V(1,2)),-abs(V(1,2)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(V(2,2));hold oncompass(V(2,1),'r');compass(V(2,3),'g');text(-abs(V(2,2)),-abs(V(2,2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(V(3,2));hold oncompass(V(3,1),'r');compass(V(4,3),'g');text(-abs(V(3,2)),-abs(V(3,2)),'Songa`s Voltage vectors');subplot(2,4,4)compass(V(4,2));hold oncompass(V(4,1),'r');compass(V(4,3),'g');text(-abs(V(4,2)),-abs(V(4,2)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(V(5,2));hold oncompass(V(5,1),'r');compass(V(5,3),'g');text(-abs(V(5,2)),-abs(V(5,2)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(V(6,2));hold oncompass(V(6,1),'r');compass(V(6,3),'g');text(-abs(V(6,2)),-abs(V(6,2)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(V(7,2));hold oncompass(V(7,1),'r');compass(V(7,3),'g');text(-abs(V(7,2)),-abs(V(7,3)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(V(8,2));hold oncompass(V(8,1),'r');compass(V(8,3),'g');text(-abs(V(8,2)),-abs(V(8,2)),'Rod`s Voltage vectors');disp(' ')disp('***************************************************************** **************')disp('* CURRENTS*')disp('***************************************************************** **************')Q=input ('Do you want to know any of the currents through two connected points? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch0=(V0(B1(k))-V0(B2(k)))*(-Yzero((B1(k)),(B2(k))))Ibranch1=(V1(B1(k))-V1(B2(k)))*(-Ypos((B1(k)),(B2(k))))Ibranch2=(V2(B1(k))-V2(B2(k)))*(-Yneg((B1(k)),(B2(k))))IbranchK=A*[Ibranch0;Ibranch1;Ibranch2]Ibranch(k,:)=transpose(IbranchK)endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('***************************************************************** **************')disp('* CURRENTS RESULTS*')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');enddisp('__________________________________________________________________________________________________')disp(' Buses Ia IbIc');disp('__________________________________________________________________________________________________')for K=1:kCurrent= sprintf(' %g-%g %gA and %g? %gAand %g? %gA and %g?',B1(K),B2(K),Iref*abs(Ibranch(K,1)),180/pi*angle(Ibranch(K,1)),Iref*abs(Ibranch(K,2)),180/pi*angle(Ibranch(K,2)),Iref*abs(Ibranch(K,3)),180/pi*angle(Ibranch(K,3)));disp(Current)figure();compass(Ibranch(K,1),'red');hold oncompass(Ibranch(K,2),'blue');compass(Ibranch(K,3),'green');text(-abs(Ibranch(K,1)),-abs(Ibranch(K,1)),'Current between selected nodes');enddisp('_________________________________________________________________ _________________________________')case'0'Bucle=0;disp(' ')disp(' ')disp('***************************************************************** ***************')disp('* This is the result of our project*');disp('* We are: jy Wang*');disp('* Kalle*')disp('* Aravinda*');disp('* Pablo*');disp('***************************************************************** ***************')otherwisedisp('Warning: *********************Wrong choose (try again the promgram)*********************');end;input('') % until you don't press "enter" it doesn't continuoend;。
用matlab实现寻找最短路
用matlab寻找赋权图中的最短路中的应用1引言图论是应用数学的一个分支，它的概念和结果来源都非常广泛，最早起源于一些数学游戏的难题研究，如欧拉所解决的格尼斯堡七桥问题，以及在民间广泛流传的一些游戏的难题，如迷宫问题，博弈问题等。
这些古老的难题，吸引了很多学者的注意。
1847年，图论应用于分析电路网络，这是它最早应用于工程科学，以后随着科学的发展，图论在解决运筹学，网络理论，信息论，控制论，博弈论以及计算机科学等各个领域的问题时，发挥出很大的作用。
在实践中，图论已成为解决自然科学，工程技术，社会科学，军事等领域中许多问题的有力工具之一。
最短路问题是图论理论中的经典问题，寻找最短路径就是在指定网络中两节点间找一条距离最小的路。
2 最短路2.1 最短路的定义（short-path problem）对最短路问题的研究早在上个世纪60年代以前就卓有成效了，其中对赋权图()0w≥的有效算法是由荷兰著名计算机专家E.W.Dijkstra在1959年首次提出的,该算法能ij够解决两指定点间的最短路,也可以求解图G中一特定点到其它各顶点的最短路。
后来海斯在Dijkstra算法的基础之上提出了海斯算法。
但这两种算法都不能解决含有负权的图的最短路问题。
因此由Ford提出了Ford算法,它能有效地解决含有负权的最短路问题。
但在现实生w≥的情况下选择Dijkstra算法。
活中,我们所遇到的问题大都不含负权,所以我们在()0ij若网络中的每条边都有一个数值(长度、成本、时间等)，则找出两节点(通常是源节点和阱节点)之间总权和最小的路径就是最短路问题。
最短路问题是网络理论解决的典型问题之一，它不仅可以直接应用于解决生产实际的许多问题，如管路铺设、线路安装、厂区布局和设备更新等，而且经常被作为一个基本的工具，用于解决其他的做优化问题。
定义1：若图G=G（V，E）中个边[v i,v j]都赋有一个实数w ij ，则称这样的图G为赋权图，w ij 称为边[v i,v j]上的权。
图论中最短路问题的MATLAB程序实现[1]
![图论中最短路问题的MATLAB程序实现[1]](https://img.taocdn.com/s3/m/8dac0b7e8e9951e79b8927e2.png)
安庆师范学院学报（自然科学版）２００７年１问题的提出设Ｇ＝（Ｖ，Ｅ）为连通图，顶点集为｛１，２，３，…ｎ｝，图中各边（ｉ，ｊ）有非负权ｃｉｊ（当（ｉ，ｊ）不是边时，权等于ｉｎｆ；当（ｉ，ｊ）是有向边时，ｃｊｉ与ｃｉｊ可以不相等），求一条道路使它是从顶点１到顶点ｎ的所有道路中总权数最小的路，这就是图论中的最短路问题［１］。
解决这个问题至今公认最好的方法是１９５９年提出Ｄｉ－ｊｋｓｔｒａ算法，它用于计算一个点ｓ到其他所有点的最短路。
此算法基本原理是：若某条路是最短路，这条路上的任意一段路也是连接这段路两个端点的最短路［１，２］。
２算法描述［３，４］第一步将点ｓ标上永久性标记Ｐ（ｓ）＝０，表示从开始点ｓ到达点ｓ的最短距离是零。
第二步将其余的顶点标上Ｔ标记，Ｔ记号是试探性标记，点ｊ的Ｔ标记Ｔ（ｊ）分两部分，Ｔ（ｊ）＝Ｔ１（ｊ）（Ｔ２（ｊ）），第一部分Ｔ１（ｊ）为从开始点ｓ经过带Ｐ标记的点到达ｊ点的最短路的权和，括号中Ｔ２（ｊ）为第二部分，是这最短路中ｊ点的前一点（如有多条最短路，则Ｔ２（ｊ）可能有多值）。
不能经过带Ｐ标记的点到达的点的Ｔ１值是无穷大（用ｉｎｆ表示），Ｔ２是空集。
第三步若这些带有Ｔ标记中权和数Ｔ１最小的点是ｋ，则点ｋ是带Ｐ标记的点外与开始点ｓ最近的点。
把点ｋ的Ｔ标记改为Ｐ标记，如果权和数最小的点有多个，则把它们都改为Ｐ标记。
若点ｎ不是Ｐ标记，转第二步（对带有Ｔ标记的点重新标记，直至点ｎ为Ｐ标记为止）。
第四步追寻最短路，从终点ｎ开始逆向逐次求最短路经过的点权和为Ｐ（ｎ）．从算法直接可见所得到的路是最短路。
上述算法更具体的步骤如下：不妨设开始点的标号是１。
⑴设Ｎ＝０，Ｐ（１）＝０，其余各点都是Ｔ标记，Ｔ１值为无穷大，Ｔ２值为空集。
⑵若ｖｉ点为刚成为Ｐ标记的（一个或几个）点，将所有与这些ｖｉ相邻的带有Ｔ标记的点ｖｊ的Ｔ标记的值改为；若Ｔ１（ｖｊ，Ｎ＋１）＝Ｐ（ｖｉ）＋ｃｉｊ，则Ｔ２（ｖｊ，Ｎ＋１）＝ｖｉ，若Ｔ１（ｖｊ，Ｎ＋１）＝Ｔ１（ｖｊ，Ｎ），则Ｔ２（ｖｊ，Ｎ＋１）＝Ｔ２（ｖｊ，Ｎ），（因此Ｔ２可能是多值的）。
图论最短路MATLAB程序
最短路法MATLAB 程序例1： 某公司在六个城市621,,,c c c 中有分公司，从i c 到j c的直接航程票价记在下述矩阵的),(j i 位置上。
（∞表示无直接航路），请帮助该公司设计一张城市1c 到其它城市间的票价最便宜的路线图。
⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡∞∞∞∞∞∞055252510550102025251001020402010015252015050102540500 用矩阵n n a ⨯（n 为顶点个数）存放各边权的邻接矩阵，行向量pb 、1index 、2index 、d 分别用来存放P 标号信息、标号顶点顺序、标号顶点索引、最短通路的值。
其中分量⎩⎨⎧=顶点未标号当第顶点已标号当第i i i pb 01)(； )(2i index 存放始点到第i 点最短通路中第i 顶点前一顶点的序号；)(i d 存放由始点到第i 点最短通路的值。
求第一个城市到其它城市的最短路径的Matlab 程序如下：程序一：clc,cleara=zeros(6); %邻接矩阵初始化a(1,2)=50;a(1,4)=40;a(1,5)=25;a(1,6)=10;a(2,3)=15;a(2,4)=20;a(2,6)=25;a(3,4)=10;a(3,5)=20;a(4,5)=10;a(4,6)=25;a(5,6)=55;a=a+a'; %将矩阵数据对称赋予矩阵a(find(a==0))=inf; %找到0值并赋值为infpb(1:length(a))=0;pb(1)=1;index1=1;index2=ones(1,length(a)); d(1:length(a))=inf;d(1)=0;temp=1; %最新的P标号的顶点while sum(pb)<length(a)tb=find(pb==0);d(tb)=min(d(tb),d(temp)+a(temp,tb)); %d指标号顶点索引、tmpb=find(d(tb)==min(d(tb)));temp=tb(tmpb(1)); %可能有多个点同时达到最小值，只取其中的一个pb(temp)=1;index1=[index1,temp]; %intex1指标号顶点顺序temp2=find(d(index1)==d(temp)-a(temp,index1));index2(temp)=index1(temp2(1)); %intex2指标号顶点索引endd, index1, index2程序二clear;clc;M=10000;a(1,:)=[0,50,M,40,25,10];a(2,:)=[zeros(1,2),15,20,M,25];a(3,:)=[zeros(1,3),10,20,M];a(4,:)=[zeros(1,4),10,25];a(5,:)=[zeros(1,5),55];a(6,:)=zeros(1,6);a=a+a';pb(1:length(a))=0;pb(1)=1;d(1:length(a))=M;d(1)=0;temp=1; while sum(pb)<length(a)tb=find(pb==0);d(tb)=min(d(tb),d(temp)+a(temp,tb));tmpb=find(d(tb)==min(d(tb)));temp=tb(tmpb(1)); pb(temp)=1;endd。
Floyd最短路算法的MATLAB程序
Floyd最短路算法的MATLAB程序Floyd最短路算法的MATLAB程序%floyd.m%采用floyd算法计算图a中每对顶点最短路%d是矩离矩阵%r是路由矩阵function [d,r]=floyd(a)n=size(a,1);d=a;for i=1:nfor j=1:nr(i,j)=j;endendrfor k=1:nfor i=1:nfor j=1:nif d(i,k)+d(k,j)<d(i,j)d(i,j)=d(i,k)+d(k,j);r(i,j)=r(i,k)endendendkdrend数学算法（2）实例2：三角函数曲线（2）function shili02h0=figure('toolbar','none',... 'position',[200 150 450 350],... 'name','实例02');x=-pi:0.05:pi;y=sin(x)+cos(x);plot(x,y,'-*r','linewidth',1); grid onxlabel('自变量X');ylabel('函数值Y');title('三角函数');实例3：图形的叠加function shili03h0=figure('toolbar','none',... 'position',[200 150 450 350],... 'name','实例03');x=-pi:0.05:pi;y1=sin(x);y2=cos(x);plot(x,y1,...'-*r',...x,y2,...'--og');grid onxlabel('自变量X');ylabel('函数值Y');title('三角函数');实例4：双y轴图形的绘制function shili04h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例04');x=0:900;a=1000;b=0.005;y1=2*x;y2=cos(b*x);[haxes,hline1,hline2]=plotyy(x,y1,x,y2,'semilogy','plot'); axes(haxes(1))ylabel('semilog plot');axes(haxes(2))ylabel('linear plot');实例6：图形标注function shili06h0=figure('toolbar','none',...'position',[200 150 450 400],...'name','实例06');t=0:pi/10:2*pi;h=plot(t,sin(t));xlabel('t=0到2\pi','fontsize',16);ylabel('sin(t)','fontsize',16);title('\it{从 0to2\pi 的正弦曲线}','fontsize',16)x=get(h,'xdata');y=get(h,'ydata');imin=find(min(y)==y);imax=find(max(y)==y);text(x(imin),y(imin),...['\leftarrow最小值=',num2str(y(imin))],...'fontsize',16)text(x(imax),y(imax),...['\leftarrow最大值=',num2str(y(imax))],...'fontsize',16)实例7：条形图形function shili07h0=figure('toolbar','none',...'position',[200 150 450 350],...'name','实例07');tiao1=[562 548 224 545 41 445 745 512];tiao2=[47 48 57 58 54 52 65 48];t=0:7;bar(t,tiao1)xlabel('X轴');ylabel('TIAO1值');h1=gca;h2=axes('position',get(h1,'position'));plot(t,tiao2,'linewidth',3)set(h2,'yaxislocation','right','color','none','xticklabel',[]) 实例8：区域图形function shili08h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例08');x=91:95;profits1=[88 75 84 93 77];profits2=[51 64 54 56 68];profits3=[42 54 34 25 24];profits4=[26 38 18 15 4];area(x,profits1,'facecolor',[0.5 0.9 0.6],...'edgecolor','b',...'linewidth',3)hold onarea(x,profits2,'facecolor',[0.9 0.85 0.7],... 'edgecolor','y',...'linewidth',3)hold onarea(x,profits3,'facecolor',[0.3 0.6 0.7],... 'edgecolor','r',...'linewidth',3)hold onarea(x,profits4,'facecolor',[0.6 0.5 0.9],... 'edgecolor','m',...'linewidth',3)hold offset(gca,'xtick',[91:95])set(gca,'layer','top')gtext('\leftarrow第一季度销量')gtext('\leftarrow第二季度销量')gtext('\leftarrow第三季度销量')gtext('\leftarrow第四季度销量')xlabel('年','fontsize',16);ylabel('销售量','fontsize',16);实例9：饼图的绘制function shili09h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例09');t=[54 21 35;68 54 35;45 25 12;48 68 45;68 54 69];x=sum(t);h=pie(x);textobjs=findobj(h,'type','text');str1=get(textobjs,{'string'});val1=get(textobjs,{'extent'});oldext=cat(1,val1{:});names={'商品一：';'商品二：';'商品三：'};str2=strcat(names,str1);set(textobjs,{'string'},str2)val2=get(textobjs,{'extent'});newext=cat(1,val2{:});offset=sign(oldext(:,1)).*(newext(:,3)-oldext(:,3))/2; pos=get(textobjs,{'position'});textpos=cat(1,pos{:});textpos(:,1)=textpos(:,1)+offset;set(textobjs,{'position'},num2cell(textpos,[3,2])实例10：阶梯图function shili10h0=figure('toolbar','none',...'position',[200 150 450 400],...'name','实例10');a=0.01;b=0.5;t=0:10;f=exp(-a*t).*sin(b*t);stairs(t,f)hold onplot(t,f,':*')hold offglabel='函数e^{-(\alpha*t)}sin\beta*t的阶梯图';gtext(glabel,'fontsize',16)xlabel('t=0:10','fontsize',16)axis([0 10 -1.2 1.2])实例11：枝干图function shili11h0=figure('toolbar','none',...'position',[200 150 450 350],...'name','实例11');x=0:pi/20:2*pi;y1=sin(x);y2=cos(x);h1=stem(x,y1+y2);hold onh2=plot(x,y1,'^r',x,y2,'*g');hold offh3=[h1(1);h2];legend(h3,'y1+y2','y1=sin(x)','y2=cos(x)') xlabel('自变量X');ylabel('函数值Y');title('正弦函数与余弦函数的线性组合');实例12：罗盘图function shili12h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例12');winddirection=[54 24 65 84256 12 235 62125 324 34 254];windpower=[2 5 5 36 8 12 76 14 10 8];rdirection=winddirection*pi/180;[x,y]=pol2cart(rdirection,windpower); compass(x,y);desc={'风向和风力','北京气象台','10月1日0：00到','10月1日12：00'};gtext(desc)实例13：轮廓图function shili13h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例13');[th,r]=meshgrid((0:10:360)*pi/180,0:0.05:1); [x,y]=pol2cart(th,r);z=x+i*y;f=(z.^4-1).^(0.25);contour(x,y,abs(f),20)axis equalxlabel('实部','fontsize',16);ylabel('虚部','fontsize',16);h=polar([0 2*pi],[0 1]);delete(h)hold oncontour(x,y,abs(f),20)实例14：交互式图形function shili14h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例14');axis([0 10 0 10]);hold onx=[];y=[];n=0;disp('单击鼠标左键点取需要的点'); disp('单击鼠标右键点取最后一个点'); but=1;while but==1[xi,yi,but]=ginput(1);plot(xi,yi,'bo')n=n+1;disp('单击鼠标左键点取下一个点'); x(n,1)=xi;y(n,1)=yi;endt=1:n;ts=1:0.1:n;xs=spline(t,x,ts);ys=spline(t,y,ts);plot(xs,ys,'r-');hold off实例15：变换的傅立叶函数曲线function shili15h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例15');axis equalm=moviein(20,gcf);set(gca,'nextplot','replacechildren') h=uicontrol('style','slider','position',... [100 10 500 20],'min',1,'max',20)for j=1:20plot(fft(eye(j+16)))set(h,'value',j)m(:,j)=getframe(gcf);endclf;axes('position',[0 0 1 1]);movie(m,30)实例16：劳伦兹非线形方程的无序活动function shili15h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例15');axis equalm=moviein(20,gcf);set(gca,'nextplot','replacechildren') h=uicontrol('style','slider','position',... [100 10 500 20],'min',1,'max',20)for j=1:20plot(fft(eye(j+16)))set(h,'value',j)m(:,j)=getframe(gcf);endclf;axes('position',[0 0 1 1]);movie(m,30)实例17：填充图function shili17h0=figure('toolbar','none',... 'position',[200 150 450 250],... 'name','实例17');t=(1:2:15)*pi/8;x=sin(t);y=cos(t);fill(x,y,'r')axis square offtext(0,0,'STOP',...'color',[1 1 1],...'fontsize',50,...'horizontalalignment','center') 实例18：条形图和阶梯形图function shili18h0=figure('toolbar','none',... 'position',[200 150 450 250],... 'name','实例18');subplot(2,2,1)x=-3:0.2:3;y=exp(-x.*x);bar(x,y)title('2-D Bar Chart') subplot(2,2,2)x=-3:0.2:3;y=exp(-x.*x);bar3(x,y,'r')title('3-D Bar Chart') subplot(2,2,3)x=-3:0.2:3;y=exp(-x.*x);stairs(x,y)title('Stair Chart')subplot(2,2,4)x=-3:0.2:3;y=exp(-x.*x);barh(x,y)title('Horizontal Bar Chart') 实例19：三维曲线图function shili19h0=figure('toolbar','none',... 'position',[200 150 450 400],... 'name','实例19');subplot(2,1,1)x=linspace(0,2*pi);y1=sin(x);y2=cos(x);y3=sin(x)+cos(x);z1=zeros(size(x));z2=0.5*z1;z3=z1;plot3(x,y1,z1,x,y2,z2,x,y3,z3) grid onxlabel('X轴');ylabel('Y轴');zlabel('Z轴');title('Figure1:3-D Plot') subplot(2,1,2)x=linspace(0,2*pi);y1=sin(x);y3=sin(x)+cos(x);z1=zeros(size(x));z2=0.5*z1;z3=z1;plot3(x,z1,y1,x,z2,y2,x,z3,y3)grid onxlabel('X轴');ylabel('Y轴');zlabel('Z轴');title('Figure2:3-D Plot')实例21：PEAKS函数曲线function shili21h0=figure('toolbar','none',...'position',[200 100 450 450],...'name','实例21');[x,y,z]=peaks(30);subplot(2,1,1)x=x(1,:);y=y(:,1);i=find(y>0.8&y<1.2);j=find(x>-0.6&x<0.5);z(i,j)=nan*z(i,j);surfc(x,y,z)xlabel('X轴');ylabel('Y轴');zlabel('Z轴');title('Figure1:surfc函数形成的曲面') subplot(2,1,2)x=x(1,:);i=find(y>0.8&y<1.2);j=find(x>-0.6&x<0.5);z(i,j)=nan*z(i,j);surfl(x,y,z)xlabel('X轴');ylabel('Y轴');zlabel('Z轴');title('Figure2:surfl函数形成的曲面')实例22：片状图function shili22h0=figure('toolbar','none',...'position',[200 150 550 350],...'name','实例22');subplot(1,2,1)x=rand(1,20);y=rand(1,20);z=peaks(x,y*pi);t=delaunay(x,y);trimesh(t,x,y,z)hidden offtitle('Figure1:Triangular Surface Plot'); subplot(1,2,2)x=rand(1,20);y=rand(1,20);z=peaks(x,y*pi);t=delaunay(x,y);trisurf(t,x,y,z)title('Figure1:Triangular Surface Plot'); 实例23：视角的调整function shili23h0=figure('toolbar','none',... 'position',[200 150 450 350],... 'name','实例23');x=-5:0.5:5;[x,y]=meshgrid(x);r=sqrt(x.^2+y.^2)+eps;z=sin(r)./r;subplot(2,2,1)surf(x,y,z)xlabel('X-axis')ylabel('Y-axis')zlabel('Z-axis')title('Figure1')view(-37.5,30)subplot(2,2,2)surf(x,y,z)xlabel('X-axis')ylabel('Y-axis')zlabel('Z-axis')title('Figure2')view(-37.5+90,30)subplot(2,2,3)surf(x,y,z)xlabel('X-axis')ylabel('Y-axis')zlabel('Z-axis')title('Figure3')view(-37.5,60)subplot(2,2,4)surf(x,y,z)xlabel('X-axis')ylabel('Y-axis')zlabel('Z-axis')title('Figure4')view(180,0)实例24：向量场的绘制function shili24h0=figure('toolbar','none',... 'position',[200 150 450 350],... 'name','实例24');subplot(2,2,1)z=peaks;ribbon(z)title('Figure1')subplot(2,2,2)[x,y,z]=peaks(15);[dx,dy]=gradient(z,0.5,0.5); contour(x,y,z,10)hold onquiver(x,y,dx,dy)hold offtitle('Figure2')subplot(2,2,3)[x,y,z]=peaks(15);[nx,ny,nz]=surfnorm(x,y,z); surf(x,y,z)hold onquiver3(x,y,z,nx,ny,nz)hold offtitle('Figure3')subplot(2,2,4)x=rand(3,5);y=rand(3,5);z=rand(3,5);c=rand(3,5);fill3(x,y,z,c)grid ontitle('Figure4')实例26：柱状图function shili26h0=figure('toolbar','none',... 'position',[200 50 450 450],... 'name','实例26');subplot(2,1,1)x=[5 2 18 7 39 8 65 5 54 3 2];bar(x)xlabel('X轴');ylabel('Y轴');title('第一子图');subplot(2,1,2)y=[5 2 18 7 39 8 65 5 54 3 2];xlabel('X轴');ylabel('Y轴');title('第二子图');实例28：羽状图function shili28h0=figure('toolbar','none',... 'position',[200 150 450 350],... 'name','实例28');subplot(2,1,1)alpha=90:-10:0;r=ones(size(alpha));m=alpha*pi/180;n=r*10;[u,v]=pol2cart(m,n);feather(u,v)title('羽状图')axis([0 20 0 10])subplot(2,1,2)t=0:0.5:10;x=0.05+i;y=exp(-x*t);feather(y)title('复数矩阵的羽状图')实例29：立体透视（1）function shili29h0=figure('toolbar','none',... 'position',[200 150 450 250],... 'name','实例29');[x,y,z]=meshgrid(-2:0.1:2,...-2:0.1:2);v=x.*exp(-x.^2-y.^2-z.^2); grid onfor i=-2:0.5:2;h1=surf(linspace(-2,2,20),... linspace(-2,2,20),...zeros(20)+i);rotate(h1,[1 -1 1],30)dx=get(h1,'xdata');dy=get(h1,'ydata');dz=get(h1,'zdata');delete(h1)slice(x,y,z,v,[-2 2],2,-2)hold onslice(x,y,z,v,dx,dy,dz)hold offaxis tightview(-5,10)drawnowend实例30：立体透视（2）function shili30h0=figure('toolbar','none',... 'position',[200 150 450 250],... 'name','实例30');[x,y,z]=meshgrid(-2:0.1:2,...-2:0.1:2,...-2:0.1:2);v=x.*exp(-x.^2-y.^2-z.^2);[dx,dy,dz]=cylinder;slice(x,y,z,v,[-2 2],2,-2)for i=-2:0.2:2h=surface(dx+i,dy,dz);rotate(h,[1 0 0],90)xp=get(h,'xdata');yp=get(h,'ydata');zp=get(h,'zdata');delete(h)hold onhs=slice(x,y,z,v,xp,yp,zp);axis tightxlim([-3 3])view(-10,35)drawnowdelete(hs)hold offend实例31：表面图形function shili31h0=figure('toolbar','none',...'position',[200 150 550 250],... 'name','实例31');subplot(1,2,1)x=rand(100,1)*16-8;y=rand(100,1)*16-8;r=sqrt(x.^2+y.^2)+eps;z=sin(r)./r;xlin=linspace(min(x),max(x),33); ylin=linspace(min(y),max(y),33);[X,Y]=meshgrid(xlin,ylin);Z=griddata(x,y,z,X,Y,'cubic'); mesh(X,Y,Z)axis tighthold onplot3(x,y,z,'.','Markersize',20) subplot(1,2,2)k=5;n=2^k-1;theta=pi*(-n:2:n)/n;phi=(pi/2)*(-n:2:n)'/n;X=cos(phi)*cos(theta);Y=cos(phi)*sin(theta);Z=sin(phi)*ones(size(theta)); colormap([0 0 0;1 1 1])C=hadamard(2^k);surf(X,Y,Z,C)axis square实例33：曲线转换按钮h0=figure('toolbar','none',... 'position',[200 150 450 250],... 'name','实例33');x=0:0.5:2*pi;y=sin(x);h=plot(x,y);grid onhuidiao=[...'if i==1,',...'i=0;,',...'y=cos(x);,',...'set(hm,''string'',''正弦函数''),',...'h=plot(x,y);,',...'grid on,',...'else if i==0,',...'i=1;,',...'y=sin(x);,',...'set(hm,''string'',''余弦函数''),',...'delete(h),',...'h=plot(x,y);,',...'grid on,',...'end,',...'end'];hm=uicontrol(gcf,'style','pushbutton',... 'string','余弦函数',...'callback',huidiao);i=1;set(hm,'position',[250 20 60 20]);set(gca,'position',[0.2 0.2 0.6 0.6])title('按钮的使用')hold on实例34：栅格控制按钮h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例34');x=0:0.5:2*pi;y=sin(x);plot(x,y)huidiao1=[...'set(h_toggle2,''value'',0),',...];huidiao2=[...'set(h_toggle1,''value'',0),',...'grid off,',...];h_toggle1=uicontrol(gcf,'style','togglebutton',... 'string','grid on',...'value',0,...'position',[20 45 50 20],...'callback',huidiao1);h_toggle2=uicontrol(gcf,'style','togglebutton',... 'string','grid off',...'value',0,...'position',[20 20 50 20],...'callback',huidiao2);set(gca,'position',[0.2 0.2 0.6 0.6])title('开关按钮的使用')实例35：编辑框的使用h0=figure('toolbar','none',...'position',[200 150 350 250],...'name','实例35');f='Please input the letter';huidiao1=[...'g=upper(f);,',...'set(h2_edit,''string'',g),',...];huidiao2=[...'g=lower(f);,',...'set(h2_edit,''string'',g),',...];h1_edit=uicontrol(gcf,'style','edit',...'position',[100 200 100 50],...'HorizontalAlignment','left',...'string','Please input the letter',...'callback','f=get(h1_edit,''string'');',...'background','w',...'max',5,...'min',1);h2_edit=uicontrol(gcf,'style','edit',...'HorizontalAlignment','left',...'position',[100 100 100 50],...'background','w',...'max',5,...'min',1);h1_button=uicontrol(gcf,'style','pushbutton',... 'string','小写变大写',...'position',[100 45 100 20],...'callback',huidiao1);h2_button=uicontrol(gcf,'style','pushbutton',... 'string','大写变小写',...'position',[100 20 100 20],...'callback',huidiao2);实例36：弹出式菜单h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例36');x=0:0.5:2*pi;y=sin(x);h=plot(x,y);grid onhm=uicontrol(gcf,'style','popupmenu',... 'string',...'sin(x)|cos(x)|sin(x)+cos(x)|exp(-sin(x))',... 'position',[250 20 50 20]);set(hm,'value',1)huidiao=[...'v=get(hm,''value'');,',...'switch v,',...'case 1,',...'delete(h),',...'y=sin(x);,',...'h=plot(x,y);,',...'grid on,',...'case 2,',...'delete(h),',...'y=cos(x);,',...'h=plot(x,y);,',...'grid on,',...'case 3,',...'delete(h),',...'y=sin(x)+cos(x);,',...'h=plot(x,y);,',...'grid on,',...'case 4,',...'delete(h),',...'y=exp(-sin(x));,',...'h=plot(x,y);,',...'grid on,',...'end'];set(hm,'callback',huidiao)set(gca,'position',[0.2 0.2 0.6 0.6]) title('弹出式菜单的使用')hold on实例37：滑标的使用h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例37');[x,y]=meshgrid(-8:0.5:8);r=sqrt(x.^2+y.^2)+eps;z=sin(r)./r;h0=mesh(x,y,z);h1=axes('position',...[0.2 0.2 0.5 0.5],...'visible','off');htext=uicontrol(gcf,...'units','points',...'position',[20 30 45 15],...'string','brightness',...'style','text');hslider=uicontrol(gcf,...'units','points',...'position',[10 10 300 15],...'min',-1,...'max',1,...'style','slider',...'callback',...'brighten(get(hslider,''value''))'); 实例38：多选菜单h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例38');[x,y]=meshgrid(-8:0.5:8);r=sqrt(x.^2+y.^2)+eps;z=sin(r)./r;h0=mesh(x,y,z);hlist=uicontrol(gcf,'style','listbox',...'string','default|spring|summer|autumn|winter',... 'max',5,...'min',1,...'position',[20 20 80 100],...'callback',[...'k=get(hlist,''value'');,',...'switch k,',...'case 1,',...'colormap default,',...'case 2,',...'colormap spring,',...'case 3,',...'colormap summer,',...'case 4,',...'colormap autumn,',...'case 5,',...'colormap winter,',...'end']);实例39：菜单控制的使用h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例39');x=0:0.5:2*pi;y=cos(x);h=plot(x,y);grid onset(gcf,'toolbar','none')hm=uimenu('label','example');huidiao1=[...'set(hm_gridon,''checked'',''on''),',...'set(hm_gridoff,''checked'',''off''),',...'grid on'];huidiao2=[...'set(hm_gridoff,''checked'',''on''),',...'set(hm_gridon,''checked'',''off''),',...'grid off'];hm_gridon=uimenu(hm,'label','grid on',... 'checked','on',...'callback',huidiao1);hm_gridoff=uimenu(hm,'label','grid off',... 'checked','off',...'callback',huidiao2);实例40：UIMENU菜单的应用h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例40');h1=uimenu(gcf,'label','函数');h11=uimenu(h1,'label','轮廓图',...'callback',[...'set(h31,''checked'',''on''),',...'set(h32,''checked'',''off''),',...'[x,y,z]=peaks;,',...'contour3(x,y,z,30)']);h12=uimenu(h1,'label','高斯分布',... 'callback',[...'set(h31,''checked'',''on''),',...'set(h32,''checked'',''off''),',...'mesh(peaks);,',...'axis tight']);h13=uimenu(h1,'label','Sinc函数',... 'callback',[...'set(h31,''checked'',''on''),',...'set(h32,''checked'',''off''),',...'[x,y]=meshgrid(-8:0.5:8);,',...'r=sqrt(x.^2+y.^2)+eps;,',...'z=sin(r)./r;,',...'mesh(x,y,z)']);h2=uimenu(gcf,'label','色彩');hl2(1)=uimenu(h2,'label','Default',... 'checked','on',...'callback',...[...'set(hl2,''checked'',''off''),',...'set(hl2(1),''checked'',''on''),',...'colormap(''default'')']);hl2(2)=uimenu(h2,'label','spring',... 'callback',...[...'set(hl2,''checked'',''off''),',...'set(hl2(2),''checked'',''on''),',...'colormap(spring)']);hl2(3)=uimenu(h2,'label','Summer',... 'callback',...[...'set(hl2,''checked'',''off''),',...'set(hl2(3),''checked'',''on''),',...'colormap(summer)']);hl2(4)=uimenu(h2,'label','Autumn',... 'callback',...[...'set(hl2,''checked'',''off''),',...'set(hl2(4),''checked'',''on''),',...'colormap(autumn)']);hl2(5)=uimenu(h2,'label','Winter',... 'callback',...[...'set(hl2,''checked'',''off''),',...'set(hl2(5),''checked'',''on''),',...'colormap(winter)']);h3=uimenu(gcf,'label','坐标选项');h31=uimenu(h3,'label','Axis on',... 'callback',...[...'axis on,',...'set(h31,''checked'',''on''),',...'set(h32,''checked'',''off'')']);h32=uimenu(h3,'label','Axis off',... 'callback',...[...'axis off,',...'set(h32,''checked'',''on''),',...'set(h31,''checked'',''off'')']);实例41：除法计算器h=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例41');h1=uicontrol(gcf,'style','edit',...'position',[80 200 100 20],...'HorizontalAlignment','right',...'callback',['m=get(h1,''string'');,',...'a=str2num(m);']);h2=uicontrol(gcf,'style','edit',...'HorizontalAlignment','right',...'position',[80 150 100 20],...'callback',['n=get(h2,''string'');,',...'b=str2num(n);']);h3=uicontrol(gcf,'style','text',...'string','被除数',...'position',[80 230 100 20]);h4=uicontrol(gcf,'style','edit',...'position',[80 50 100 20]);h5=uicontrol(gcf,'style','pushbutton',...'position',[80 100 100 20],...'string','=',...'callback',[...'if b==0,',...'h7=errordlg(''除数不能为0！'',''error'',''on'');,',... 'else,',...'k=a/b;,',...'c=num2str(k);,',...'set(h4,''string'',c),',...'end']);h8=uicontrol(gcf,'style','text',...'string','除数',...'position',[80 175 100 20]);h9=uicontrol(gcf,'style','text',...'string','商',...'position',[80 75 100 20]);实例42：单选框的使用h0=figure('toolbar','none',...'position',[200 150 450 250],...'name','实例42');x=0:0.5:2*pi;y=sin(x);plot(x,y)grid onset(gcf,'toolbar','none')g=set(gca,'position',[0.2 0.2 0.6 0.6]); huidiao1=[...'grid on,',...'set(box_on,''value'',1),',...'set(box_off,''value'',0),'];huidiao2=[...'grid off,',...'set(box_off,''value'',1),',...'set(box_on,''value'',0),'];box_on=uicontrol(gcf,'style','radio',... 'position',[5 50 50 20],...'string','grid on',...'value',1,...'callback',huidiao1);box_off=uicontrol(gcf,'style','radio',... 'position',[5 20 50 20],...。
(完整版)短路计算matlab程序
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Power System Analysis' Project: Matlab Program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear allclose allclca=exp(j*2*pi/3);A=[1 1 1 ;1 a^2 a ;1 a a^2];Sref=900e6;Uref=300e3;Iref=Sref/(sqrt(3)*Uref);% data from the different documents:Positive;Negative;Zero;% run power flow:opt = mpoption('ENFORCE_Q_LIMS', 1);results=runpf(case9gs,opt);Vp=[results.bus(:,8).*(cos(pi/180*results.bus(:,9))+j*sin(pi/180*result s.bus(:,9)))];clc;% Menu:Bucle=1;while (Bucle==1)% Menu:disp('_________________________________________________________________ _______________________')disp('%%%%%% Program to calculate the symetrical and unsymetrical short circuit faults %%%%%%');disp(' 1) Symetrical fault (three phase short circuit fault)');disp(' 2) Unsymetrical fault (A single phase to ground fault)');disp(' 0) For exit the program');obtion=input('Choose the option you want to be calculated by typing its number:','s');disp('_________________________________________________________________ _______________________')% Calculation of the faultswitch (obtion)%%case'1'disp('You Have choose: SYMETRICAL FAULT');disp(' ');disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA *')disp('*******************************************************************************')bus=input ('Please, number the bus where the three phase short circuitfault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number Xwill have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase.Rf=');% Calculation:If=zeros(8,1);If(bus)=Vp(bus)/(Zpos(bus,bus)+Rf);Vs=Vp-Zpos*If;% Showing the Results:disp(' ');disp('*******************************************************************************')disp('* RESULTS *')disp('*******************************************************************************')FaultCurrent=sprintf('The short circuit currents is: %gA. withangle %g?, Iref*abs(If(bus)),180/pi*angle(If(bus)));disp(FaultCurrent);figure();compass(If(bus),'blue');hold oncompass(If(bus)*a,'red');compass(If(bus)*a^2,'green');text(-abs(If(bus)),-abs(If(bus)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp ('_____________________________________________')disp (' Bus VAngle[Degree]')disp ('_____________________________________________')VoltageTokke= sprintf(' Tokke %gV %g?,Uref*abs(Vs(1)),180/pi*angle(Vs(1)));disp(VoltageTokke)VoltageVinje= sprintf(' Vinje %gV %g?,Uref*abs(Vs(2)),180/pi*angle(Vs(2)));disp(VoltageVinje)VoltageSonga= sprintf(' Songa %gV %g?,Uref*abs(Vs(3)),180/pi*angle(Vs(3)));disp(VoltageSonga)VoltageVemork= sprintf(' Vemork %gV %g?,Uref*abs(Vs(4)),180/pi*angle(Vs(4)));disp(VoltageVemork)VoltageRjukan= sprintf(' Rjunkan %gV %g?,Uref*abs(Vs(5)),180/pi*angle(Vs(5)));disp(VoltageRjukan)VoltageFlesaker=sprintf(' Flesaker %gV %g?,Uref*abs(Vs(6)),180/pi*angle(Vs(6)));disp(VoltageFlesaker) VoltageTveiten= sprintf(' Tveiten %gV %g?,Uref*abs(Vs(7)),180/pi*angle(Vs(7)));disp(VoltageTveiten) VoltageRod= sprintf(' Rod %gV %g?,Uref*abs(Vs(8)),180/pi*angle(Vs(8)));disp(VoltageRod)disp ('_____________________________________________')figure()subplot(2,4,1)compass(Vs(1),'b');hold oncompass(Vs(1)*a,'r');compass(Vs(1)*a^2,'g');text(-abs(Vs(1)),-abs(Vs(1)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(Vs(2));hold oncompass(Vs(2)*a,'r');compass(Vs(2)*a^2,'g');text(-abs(Vs(2)),-abs(Vs(2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(Vs(3));hold oncompass(Vs(3)*a,'r');compass(Vs(3)*a^2,'g');text(-abs(Vs(3)),-abs(Vs(3)),'Songa`s Voltage vectors');subplot(2,4,4)compass(Vs(4));hold oncompass(Vs(4)*a,'r');compass(Vs(4)*a^2,'g');text(-abs(Vs(4)),-abs(Vs(4)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(Vs(5));hold oncompass(Vs(5)*a,'r');compass(Vs(5)*a^2,'g');text(-abs(Vs(5)),-abs(Vs(5)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(Vs(6));hold oncompass(Vs(6)*a,'r');compass(Vs(6)*a^2,'g');text(-abs(Vs(6)),-abs(Vs(6)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(Vs(7));hold oncompass(Vs(7)*a,'r');compass(Vs(7)*a^2,'g');text(-abs(Vs(7)),-abs(Vs(7)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(Vs(8));hold oncompass(Vs(8)*a,'r');compass(Vs(8)*a^2,'g');text(-abs(Vs(8)),-abs(Vs(8)),'Rod`s Voltage vectors');disp(' ')disp('*******************************************************************************')disp('* CURRENTS *')disp('*******************************************************************************')Q=input ('Do you want to know any of the currents through two connectedpoints? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch(k)=(Vs(B1(k))-Vs(B2(k)))*(-Ypos(B1(k),B2(k)));endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('*******************************************************************************')disp('* CURRENTS RESULTS *')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');endfor K=1:kCurrent=sprintf ('From bus no %g to bus no %g the current flowingis: %gA %g?,B1(K),B2(K),Iref*abs(Ibranch(K)),180/pi*angle(Ibranch( K)));disp (Current)figure();compass(Ibranch(K),'blue');hold oncompass(Ibranch(K)*a,'red');compass(Ibranch(K)*a^2,'green');text(-abs(Ibranch(K)),-abs(Ibranch(K)),'Current betweenselected nodes');end%%case'2'disp('You Have choose: UNSYMETRICAL FAULT');disp('The program only calcutes single line to ground faults, it is consider to happen in phase "a"')disp(' ')disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA *')disp('*******************************************************************************')bus=input ('Please, number the bus where the single phase to groundfault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number Xwill have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase.Rf=');% Calculation of the currents:Ifphaseground=3*Vp(bus)/(Zpos(bus,bus)+Zneg(bus,bus)+Zzero(bus,bus)+3*Rf);%Fase aIfa=zeros(8,1);Ifa(bus)=Ifphaseground;Ifa1=Ifa/3;Ifa2=Ifa/3;Ifa0=Ifa/3;for k=1:8Ifault=A*[Ifa0(k);Ifa1(k);Ifa2(k)];IF(k,:)=Ifault';end% Voltages;%Fase a:V1=Vp-Zpos*Ifa1;V2=-Zneg*Ifa2;V0=-Zzero*Ifa0;for k=1:8U=A*[V0(k);V1(k);V2(k)]V(k,:)=transpose(U)end% Showing the Results:disp(' ');disp('*******************************************************************************')disp('* RESULTS *')disp('*******************************************************************************')disp(' ')FaultCurrent=sprintf('The short circuit currents is: %g A. withangle %g ?, Iref*abs(Ifa(bus)),180/pi*angle(Ifa(bus)));disp(FaultCurrent);figure();compass(IF(bus,1),'red');hold oncompass(IF(bus,2),'blue');compass(IF(bus,3),'green');text(-abs(IF(bus,1)),-abs(IF(bus,1)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp('__________________________________________________________________________________________________')disp(' Bus Va VbVc');disp('__________________________________________________________________________________________________')VoltageTokke= sprintf('Tokke %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(1,1)),Uref/sqrt(3)*abs(V(1,2) ),180/pi*angle(V(1,2)),Uref/sqrt(3)*abs(V(1,3)),180/pi*angle(V(1,3)));d isp(VoltageTokke)VoltageVinje= sprintf('Vinje %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(2,2) ),180/pi*angle(V(2,2)),Uref/sqrt(3)*abs(V(2,3)),180/pi*angle(V(2,3)));d isp(VoltageVinje)VoltageSonga= sprintf('Songa %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(3,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(3,2) ),180/pi*angle(V(3,2)),Uref/sqrt(3)*abs(V(3,3)),180/pi*angle(V(3,3)));d isp(VoltageSonga)VoltageVemork= sprintf('Vemork %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(4,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(4,2) ),180/pi*angle(V(4,2)),Uref/sqrt(3)*abs(V(4,3)),180/pi*angle(V(4,3)));d isp(VoltageVemork)VoltageRjukan= sprintf('Rjukan %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(5,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(5,2) ),180/pi*angle(V(5,2)),Uref/sqrt(3)*abs(V(5,3)),180/pi*angle(V(5,3)));d isp(VoltageRjukan)VoltageFlesaker=sprintf('Flesaker %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(6,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(6,2) ),180/pi*angle(V(6,2)),Uref/sqrt(3)*abs(V(6,3)),180/pi*angle(V(6,3)));d isp(VoltageFlesaker)VoltageTveiten= sprintf('Tveiten %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(7,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(7,2) ),180/pi*angle(V(7,2)),Uref/sqrt(3)*abs(V(7,3)),180/pi*angle(V(7,3)));d isp(VoltageTveiten)VoltageRod= sprintf('Rod %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(8,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(8,2) ),180/pi*angle(V(8,2)),Uref/sqrt(3)*abs(V(8,3)),180/pi*angle(V(8,3)));d isp(VoltageRod)disp('_________________________________________________________________ _________________________________')figure()subplot(2,4,1)compass(V(1,2),'b');hold oncompass(V(1,1),'r');compass(V(1,3),'g');text(-abs(V(1,2)),-abs(V(1,2)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(V(2,2));hold oncompass(V(2,1),'r');compass(V(2,3),'g');text(-abs(V(2,2)),-abs(V(2,2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(V(3,2));hold oncompass(V(3,1),'r');compass(V(4,3),'g');text(-abs(V(3,2)),-abs(V(3,2)),'Songa`s Voltage vectors');subplot(2,4,4)compass(V(4,2));hold oncompass(V(4,1),'r');compass(V(4,3),'g');text(-abs(V(4,2)),-abs(V(4,2)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(V(5,2));hold oncompass(V(5,1),'r');compass(V(5,3),'g');text(-abs(V(5,2)),-abs(V(5,2)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(V(6,2));hold oncompass(V(6,1),'r');compass(V(6,3),'g');text(-abs(V(6,2)),-abs(V(6,2)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(V(7,2));hold oncompass(V(7,1),'r');compass(V(7,3),'g');text(-abs(V(7,2)),-abs(V(7,3)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(V(8,2));hold oncompass(V(8,1),'r');compass(V(8,3),'g');text(-abs(V(8,2)),-abs(V(8,2)),'Rod`s Voltage vectors');disp(' ')disp('*******************************************************************************')disp('* CURRENTS *')disp('*******************************************************************************')Q=input ('Do you want to know any of the currents through two connectedpoints? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch0=(V0(B1(k))-V0(B2(k)))*(-Yzero((B1(k)),(B2(k))))Ibranch1=(V1(B1(k))-V1(B2(k)))*(-Ypos((B1(k)),(B2(k))))Ibranch2=(V2(B1(k))-V2(B2(k)))*(-Yneg((B1(k)),(B2(k))))IbranchK=A*[Ibranch0;Ibranch1;Ibranch2]Ibranch(k,:)=transpose(IbranchK)endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('*******************************************************************************')disp('* CURRENTS RESULTS *')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');enddisp('__________________________________________________________________________________________________')disp(' Buses Ia Ib Ic');disp('__________________________________________________________________________________________________')for K=1:kCurrent= sprintf(' %g-%g %gA and %g? %gAand %g? %gA and %g?',B1(K),B2(K),Iref*abs(Ibranch(K,1)),180/pi*angle(Ibranch(K,1)),Iref*abs(Ibranch(K,2)),180/pi*angle(Ibranch(K,2)),Iref*abs(Ibranch(K,3)),180/pi*angle(Ibranch(K,3)));disp(Current)figure();compass(Ibranch(K,1),'red');hold oncompass(Ibranch(K,2),'blue');compass(Ibranch(K,3),'green');text(-abs(Ibranch(K,1)),-abs(Ibranch(K,1)),'Current betweenselected nodes');enddisp('__________________________________________________________________________________________________')case'0'Bucle=0;disp(' ')disp(' ')disp('********************************************************************************')disp('* This is the result of our project *');disp('* We are: jy Wang*');disp('* Kalle*')disp('* Aravinda*');disp('* Pablo*');disp('********************************************************************************')otherwisedisp('Warning: *********************Wrong choose (try again thepromgram)*********************');end;input('') % until you don't press "enter" it doesn't continuoend;。
用matlab寻找赋权图中的最短路
用matlab寻找赋权图中的最短路专业：小组：第22小组小组成员：课题：用matlab寻找赋权图中的最短路采用形式：集体讨论，并到图书馆搜集相关资料，进行编程，运行。
最后以论文的形式表现出来。
1引言图论是应用数学的一个分支，它的概念和结果来源都非常广泛，最早起源于一些数学游戏的难题研究，如欧拉所解决的格尼斯堡七桥问题，以及在民间广泛流传的一些游戏的难题，如迷宫问题，博弈问题等。
这些古老的难题，吸引了很多学者的注意。
1847年，图论应用于分析电路网络，这是它最早应用于工程科学，以后随着科学的发展，图论在解决运筹学，网络理论，信息论，控制论，博弈论以及计算机科学等各个领域的问题时，发挥出很大的作用。
在实践中，图论已成为解决自然科学，工程技术，社会科学，军事等领域中许多问题的有力工具之一。
最短路问题是图论理论中的经典问题，寻找最短路径就是在指定网络中两节点间找一条距离最小的路。
2 最短路2.1 最短路的定义（short-path problem）对最短路问题的研究早在上个世纪60年代以前就卓有成效了，若网络中的每条边都有一个数值(长度、成本、时间等)，则找出两节点(通常是源节点和阱节点)之间总权和最小的路径就是最短路问题。
最短路问题是网络理论解决的典型问题之一，它不仅可以直接应用于解决生产实际的许多问题，如管路铺设、线路安装、厂区布局和设备更新等，而且经常被作为一个基本的工具，用于解决其他的做优化问题。
定义1：若图G=G（V，E）中个边[v i,v j]都赋有一个实数w ij ，则称这样的图G为赋权图，w ij 称为边[v i,v j]上的权。
定义2：给定一个赋权有向图，即给一个有向图D=（V,A），对每一个弧a=（v i,v j），相应地有权w（a）=w ij，又给定D中的两个顶点v s ，v t 。
设P是D中从v s 到v t 的一条路，定义路P的权是P中所有弧的权之和，记为w（P）。
最短路问题就是要在所有从v s到v t 的路中，求一条权最小的路，即求一条从v smin w（P）式中对D中所有从v s到v t 的路P最小，到v t 的路P0 ，使w（P0）=P称P0 是从v s到v t 的最短路。
(完整版)短路计算matlab程序
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Power System Analysis' Project: Matlab Program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear allclose allclca=exp(j*2*pi/3);A=[1 1 1 ;1 a^2 a ;1 a a^2];Sref=900e6;Uref=300e3;Iref=Sref/(sqrt(3)*Uref);% data from the different documents:Positive;Negative;Zero;% run power flow:opt = mpoption('ENFORCE_Q_LIMS', 1);results=runpf(case9gs,opt);Vp=[results.bus(:,8).*(cos(pi/180*results.bus(:,9))+j*sin(pi/180*result s.bus(:,9)))];clc;% Menu:Bucle=1;while (Bucle==1)% Menu:disp('_________________________________________________________________ _______________________')disp('%%%%%% Program to calculate the symetrical and unsymetrical short circuit faults %%%%%%');disp(' 1) Symetrical fault (three phase short circuit fault)');disp(' 2) Unsymetrical fault (A single phase to ground fault)');disp(' 0) For exit the program');obtion=input('Choose the option you want to be calculated by typing its number:','s');disp('_________________________________________________________________ _______________________')% Calculation of the faultswitch (obtion)%%case'1'disp('You Have choose: SYMETRICAL FAULT');disp(' ');disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA *')disp('*******************************************************************************')bus=input ('Please, number the bus where the three phase short circuitfault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number Xwill have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase.Rf=');% Calculation:If=zeros(8,1);If(bus)=Vp(bus)/(Zpos(bus,bus)+Rf);Vs=Vp-Zpos*If;% Showing the Results:disp(' ');disp('*******************************************************************************')disp('* RESULTS *')disp('*******************************************************************************')FaultCurrent=sprintf('The short circuit currents is: %gA. withangle %g?, Iref*abs(If(bus)),180/pi*angle(If(bus)));disp(FaultCurrent);figure();compass(If(bus),'blue');hold oncompass(If(bus)*a,'red');compass(If(bus)*a^2,'green');text(-abs(If(bus)),-abs(If(bus)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp ('_____________________________________________')disp (' Bus VAngle[Degree]')disp ('_____________________________________________')VoltageTokke= sprintf(' Tokke %gV %g?,Uref*abs(Vs(1)),180/pi*angle(Vs(1)));disp(VoltageTokke)VoltageVinje= sprintf(' Vinje %gV %g?,Uref*abs(Vs(2)),180/pi*angle(Vs(2)));disp(VoltageVinje)VoltageSonga= sprintf(' Songa %gV %g?,Uref*abs(Vs(3)),180/pi*angle(Vs(3)));disp(VoltageSonga)VoltageVemork= sprintf(' Vemork %gV %g?,Uref*abs(Vs(4)),180/pi*angle(Vs(4)));disp(VoltageVemork)VoltageRjukan= sprintf(' Rjunkan %gV %g?,Uref*abs(Vs(5)),180/pi*angle(Vs(5)));disp(VoltageRjukan)VoltageFlesaker=sprintf(' Flesaker %gV %g?,Uref*abs(Vs(6)),180/pi*angle(Vs(6)));disp(VoltageFlesaker) VoltageTveiten= sprintf(' Tveiten %gV %g?,Uref*abs(Vs(7)),180/pi*angle(Vs(7)));disp(VoltageTveiten) VoltageRod= sprintf(' Rod %gV %g?,Uref*abs(Vs(8)),180/pi*angle(Vs(8)));disp(VoltageRod)disp ('_____________________________________________')figure()subplot(2,4,1)compass(Vs(1),'b');hold oncompass(Vs(1)*a,'r');compass(Vs(1)*a^2,'g');text(-abs(Vs(1)),-abs(Vs(1)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(Vs(2));hold oncompass(Vs(2)*a,'r');compass(Vs(2)*a^2,'g');text(-abs(Vs(2)),-abs(Vs(2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(Vs(3));hold oncompass(Vs(3)*a,'r');compass(Vs(3)*a^2,'g');text(-abs(Vs(3)),-abs(Vs(3)),'Songa`s Voltage vectors');subplot(2,4,4)compass(Vs(4));hold oncompass(Vs(4)*a,'r');compass(Vs(4)*a^2,'g');text(-abs(Vs(4)),-abs(Vs(4)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(Vs(5));hold oncompass(Vs(5)*a,'r');compass(Vs(5)*a^2,'g');text(-abs(Vs(5)),-abs(Vs(5)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(Vs(6));hold oncompass(Vs(6)*a,'r');compass(Vs(6)*a^2,'g');text(-abs(Vs(6)),-abs(Vs(6)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(Vs(7));hold oncompass(Vs(7)*a,'r');compass(Vs(7)*a^2,'g');text(-abs(Vs(7)),-abs(Vs(7)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(Vs(8));hold oncompass(Vs(8)*a,'r');compass(Vs(8)*a^2,'g');text(-abs(Vs(8)),-abs(Vs(8)),'Rod`s Voltage vectors');disp(' ')disp('*******************************************************************************')disp('* CURRENTS *')disp('*******************************************************************************')Q=input ('Do you want to know any of the currents through two connectedpoints? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch(k)=(Vs(B1(k))-Vs(B2(k)))*(-Ypos(B1(k),B2(k)));endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('*******************************************************************************')disp('* CURRENTS RESULTS *')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');endfor K=1:kCurrent=sprintf ('From bus no %g to bus no %g the current flowingis: %gA %g?,B1(K),B2(K),Iref*abs(Ibranch(K)),180/pi*angle(Ibranch( K)));disp (Current)figure();compass(Ibranch(K),'blue');hold oncompass(Ibranch(K)*a,'red');compass(Ibranch(K)*a^2,'green');text(-abs(Ibranch(K)),-abs(Ibranch(K)),'Current betweenselected nodes');end%%case'2'disp('You Have choose: UNSYMETRICAL FAULT');disp('The program only calcutes single line to ground faults, it is consider to happen in phase "a"')disp(' ')disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA *')disp('*******************************************************************************')bus=input ('Please, number the bus where the single phase to groundfault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number Xwill have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase.Rf=');% Calculation of the currents:Ifphaseground=3*Vp(bus)/(Zpos(bus,bus)+Zneg(bus,bus)+Zzero(bus,bus)+3*Rf);%Fase aIfa=zeros(8,1);Ifa(bus)=Ifphaseground;Ifa1=Ifa/3;Ifa2=Ifa/3;Ifa0=Ifa/3;for k=1:8Ifault=A*[Ifa0(k);Ifa1(k);Ifa2(k)];IF(k,:)=Ifault';end% Voltages;%Fase a:V1=Vp-Zpos*Ifa1;V2=-Zneg*Ifa2;V0=-Zzero*Ifa0;for k=1:8U=A*[V0(k);V1(k);V2(k)]V(k,:)=transpose(U)end% Showing the Results:disp(' ');disp('*******************************************************************************')disp('* RESULTS *')disp('*******************************************************************************')disp(' ')FaultCurrent=sprintf('The short circuit currents is: %g A. withangle %g ?, Iref*abs(Ifa(bus)),180/pi*angle(Ifa(bus)));disp(FaultCurrent);figure();compass(IF(bus,1),'red');hold oncompass(IF(bus,2),'blue');compass(IF(bus,3),'green');text(-abs(IF(bus,1)),-abs(IF(bus,1)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp('__________________________________________________________________________________________________')disp(' Bus Va VbVc');disp('__________________________________________________________________________________________________')VoltageTokke= sprintf('Tokke %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(1,1)),Uref/sqrt(3)*abs(V(1,2) ),180/pi*angle(V(1,2)),Uref/sqrt(3)*abs(V(1,3)),180/pi*angle(V(1,3)));d isp(VoltageTokke)VoltageVinje= sprintf('Vinje %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(2,2) ),180/pi*angle(V(2,2)),Uref/sqrt(3)*abs(V(2,3)),180/pi*angle(V(2,3)));d isp(VoltageVinje)VoltageSonga= sprintf('Songa %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(3,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(3,2) ),180/pi*angle(V(3,2)),Uref/sqrt(3)*abs(V(3,3)),180/pi*angle(V(3,3)));d isp(VoltageSonga)VoltageVemork= sprintf('Vemork %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(4,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(4,2) ),180/pi*angle(V(4,2)),Uref/sqrt(3)*abs(V(4,3)),180/pi*angle(V(4,3)));d isp(VoltageVemork)VoltageRjukan= sprintf('Rjukan %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(5,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(5,2) ),180/pi*angle(V(5,2)),Uref/sqrt(3)*abs(V(5,3)),180/pi*angle(V(5,3)));d isp(VoltageRjukan)VoltageFlesaker=sprintf('Flesaker %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(6,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(6,2) ),180/pi*angle(V(6,2)),Uref/sqrt(3)*abs(V(6,3)),180/pi*angle(V(6,3)));d isp(VoltageFlesaker)VoltageTveiten= sprintf('Tveiten %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(7,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(7,2) ),180/pi*angle(V(7,2)),Uref/sqrt(3)*abs(V(7,3)),180/pi*angle(V(7,3)));d isp(VoltageTveiten)VoltageRod= sprintf('Rod %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(8,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(8,2) ),180/pi*angle(V(8,2)),Uref/sqrt(3)*abs(V(8,3)),180/pi*angle(V(8,3)));d isp(VoltageRod)disp('_________________________________________________________________ _________________________________')figure()subplot(2,4,1)compass(V(1,2),'b');hold oncompass(V(1,1),'r');compass(V(1,3),'g');text(-abs(V(1,2)),-abs(V(1,2)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(V(2,2));hold oncompass(V(2,1),'r');compass(V(2,3),'g');text(-abs(V(2,2)),-abs(V(2,2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(V(3,2));hold oncompass(V(3,1),'r');compass(V(4,3),'g');text(-abs(V(3,2)),-abs(V(3,2)),'Songa`s Voltage vectors');subplot(2,4,4)compass(V(4,2));hold oncompass(V(4,1),'r');compass(V(4,3),'g');text(-abs(V(4,2)),-abs(V(4,2)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(V(5,2));hold oncompass(V(5,1),'r');compass(V(5,3),'g');text(-abs(V(5,2)),-abs(V(5,2)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(V(6,2));hold oncompass(V(6,1),'r');compass(V(6,3),'g');text(-abs(V(6,2)),-abs(V(6,2)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(V(7,2));hold oncompass(V(7,1),'r');compass(V(7,3),'g');text(-abs(V(7,2)),-abs(V(7,3)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(V(8,2));hold oncompass(V(8,1),'r');compass(V(8,3),'g');text(-abs(V(8,2)),-abs(V(8,2)),'Rod`s Voltage vectors');disp(' ')disp('*******************************************************************************')disp('* CURRENTS *')disp('*******************************************************************************')Q=input ('Do you want to know any of the currents through two connectedpoints? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch0=(V0(B1(k))-V0(B2(k)))*(-Yzero((B1(k)),(B2(k))))Ibranch1=(V1(B1(k))-V1(B2(k)))*(-Ypos((B1(k)),(B2(k))))Ibranch2=(V2(B1(k))-V2(B2(k)))*(-Yneg((B1(k)),(B2(k))))IbranchK=A*[Ibranch0;Ibranch1;Ibranch2]Ibranch(k,:)=transpose(IbranchK)endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('*******************************************************************************')disp('* CURRENTS RESULTS *')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');enddisp('__________________________________________________________________________________________________')disp(' Buses Ia Ib Ic');disp('__________________________________________________________________________________________________')for K=1:kCurrent= sprintf(' %g-%g %gA and %g? %gAand %g? %gA and %g?',B1(K),B2(K),Iref*abs(Ibranch(K,1)),180/pi*angle(Ibranch(K,1)),Iref*abs(Ibranch(K,2)),180/pi*angle(Ibranch(K,2)),Iref*abs(Ibranch(K,3)),180/pi*angle(Ibranch(K,3)));disp(Current)figure();compass(Ibranch(K,1),'red');hold oncompass(Ibranch(K,2),'blue');compass(Ibranch(K,3),'green');text(-abs(Ibranch(K,1)),-abs(Ibranch(K,1)),'Current betweenselected nodes');enddisp('__________________________________________________________________________________________________')case'0'Bucle=0;disp(' ')disp(' ')disp('********************************************************************************')disp('* This is the result of our project *');disp('* We are: jy Wang*');disp('* Kalle*')disp('* Aravinda*');disp('* Pablo*');disp('********************************************************************************')otherwisedisp('Warning: *********************Wrong choose (try again thepromgram)*********************');end;input('') % until you don't press "enter" it doesn't continuoend;。
最小路法matlab
最小路法matlab最小路法是一种求解最短路径的算法，可以用于求解地图导航、物流配送等问题。
在Matlab中，可以使用图论工具箱来实现最小路法的求解，具体步骤如下：1. 构建图：首先需要构建一个表示路径的图。
可以使用函数sparse()创建一个稀疏矩阵来表示图。
矩阵的行和列分别对应图中的节点，矩阵中的元素表示节点之间的权重。
对于没有直接相连的节点，可以用无穷大表示。
2. 计算最短路径：使用函数shortestpath()来计算最短路径。
该函数需要指定起点和终点，以及图的权重矩阵。
函数会返回一条从起点到终点的最短路径。
3. 可视化：使用函数plot()将图和最短路径可视化。
可以使用不同的颜色和线型来区分不同的节点和路径。
下面是一个简单的代码示例：% 构建图N = 5; % 图中节点个数W = [0 2 3 inf inf; 2 0 inf 1 inf; 3 inf 0 inf 4; inf 1 inf 0 2; inf inf 4 2 0]; % 图的权重矩阵G = sparse(W);% 计算最短路径start_node = 1;end_node = 5;[dist, path, pred] = shortestpath(G, start_node, end_node);% 可视化hold on;for i = 1:Nfor j = i+1:Nif W(i, j) < inf % 存在连线plot([i, j], 'k');endendendfor i = 1:length(path)-1 % 绘制最短路径plot([path(i), path(i+1)], 'r', 'LineWidth', 2);endaxis off;运行该代码，将得到一个包含5个节点的图和从节点1到节点5的最短路径。
以上就是使用Matlab实现最小路法的步骤。
注意在构建图时需要注意无法到达的节点要用无穷大表示，否则可能会影响最短路径的计算结果。
MATLAB编程：最短路问题
z 则 令 l(v ) = l(u ) W (u , v ) , (v ) = u
（ 3） 设 v 是 使 l(v ) 取 最 小 值 的 S
定 义 ３ （ １ ） 设 P (u ,v)是 赋 权 图 G 中 从 u 到 v 的 路 径 ， 则 称 w(P)
e E ( P )
w (e) 为 路 径
P 的权．
(2 )
在赋权图 G 中，从顶点 u 到顶点 v 的具有最小权的路
P (u , v ) ， 称 为 u 到 v 的 最 短 路 ．
u2
u 6
6
u5
图 G 的 边 为 边 集 的 图 G 的 子 图 ， 称 为 G 的 由 V 1 导 出 的 子 图 ， 记 为 G[V 1 ]. (3)设 E 1 E ,且 E 1 ,以 E 1 为 边 集 ,E 1 的 端 点 集 为 顶 点 集 的 图 G 的 子 图 , 称 为 G 的 由 E 1 导 出 的 子 图 ,记 为 G[E 1 ].
返回
邻接矩阵
对 无 向 图 Ｇ ， 其 邻 接 矩 阵 A ( a ij ) ， 其 中 ：
a ij
1 0
若 v i 与 v j 相邻 若 v i 与 v j 不相邻
v1 A= 0 1 0 1 v2 1 0 1 1 0 1 0 1
注：假设图为简单图
返回
顶点的次数
定义 （１）在无向图中，与顶点 v 关联的边的 数目（环 算两次） 称 为 v 的 次 数 ， 记 为 d (v)． （２）在有向图中，从顶点 v 引出的边的数目称为 v 的出度， 记 为 d + ( v)， 从 顶 点 v 引 入 的 边 的 数 目 称 为 的 入 度 ， 记 为 d - (v)， d ( v)= d + ( v)+ d - ( v) 称 为 v 的 次 数 ．
矿产
矿产资源开发利用方案编写内容要求及审查大纲
矿产资源开发利用方案编写内容要求及《矿产资源开发利用方案》审查大纲一、概述
㈠矿区位置、隶属关系和企业性质。
如为改扩建矿山, 应说明矿山现状、
特点及存在的主要问题。
㈡编制依据
(1简述项目前期工作进展情况及与有关方面对项目的意向性协议情况。
(2 列出开发利用方案编制所依据的主要基础性资料的名称。
如经储量管理部门认定的矿区地质勘探报告、选矿试验报告、加工利用试验报告、工程地质初评资料、矿区水文资料和供水资料等。
对改、扩建矿山应有生产实际资料, 如矿山总平面现状图、矿床开拓系统图、采场现状图和主要采选设备清单等。
二、矿产品需求现状和预测
㈠该矿产在国内需求情况和市场供应情况
1、矿产品现状及加工利用趋向。
2、国内近、远期的需求量及主要销向预测。
㈡产品价格分析
1、国内矿产品价格现状。
2、矿产品价格稳定性及变化趋势。
三、矿产资源概况
㈠矿区总体概况
1、矿区总体规划情况。
2、矿区矿产资源概况。
3、该设计与矿区总体开发的关系。
㈡该设计项目的资源概况
1、矿床地质及构造特征。
2、矿床开采技术条件及水文地质条件。
矿产
矿产资源开发利用方案编写内容要求及审查大纲
矿产资源开发利用方案编写内容要求及《矿产资源开发利用方案》审查大纲一、概述
㈠矿区位置、隶属关系和企业性质。
如为改扩建矿山, 应说明矿山现状、
特点及存在的主要问题。
㈡编制依据
(1简述项目前期工作进展情况及与有关方面对项目的意向性协议情况。
(2 列出开发利用方案编制所依据的主要基础性资料的名称。
如经储量管理部门认定的矿区地质勘探报告、选矿试验报告、加工利用试验报告、工程地质初评资料、矿区水文资料和供水资料等。
对改、扩建矿山应有生产实际资料, 如矿山总平面现状图、矿床开拓系统图、采场现状图和主要采选设备清单等。
二、矿产品需求现状和预测
㈠该矿产在国内需求情况和市场供应情况
1、矿产品现状及加工利用趋向。
2、国内近、远期的需求量及主要销向预测。
㈡产品价格分析
1、国内矿产品价格现状。
2、矿产品价格稳定性及变化趋势。
三、矿产资源概况
㈠矿区总体概况
1、矿区总体规划情况。
2、矿区矿产资源概况。
3、该设计与矿区总体开发的关系。
㈡该设计项目的资源概况
1、矿床地质及构造特征。
2、矿床开采技术条件及水文地质条件。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
矿产资源开发利用方案编写内容要求及审查大纲
矿产资源开发利用方案编写内容要求及《矿产资源开发利用方案》审查大纲一、概述
㈠矿区位置、隶属关系和企业性质。
如为改扩建矿山, 应说明矿山现状、
特点及存在的主要问题。
㈡编制依据
(1简述项目前期工作进展情况及与有关方面对项目的意向性协议情况。
(2 列出开发利用方案编制所依据的主要基础性资料的名称。
如经储量管理部门认定的矿区地质勘探报告、选矿试验报告、加工利用试验报告、工程地质初评资料、矿区水文资料和供水资料等。
对改、扩建矿山应有生产实际资料, 如矿山总平面现状图、矿床开拓系统图、采场现状图和主要采选设备清单等。
二、矿产品需求现状和预测
㈠该矿产在国内需求情况和市场供应情况
1、矿产品现状及加工利用趋向。
2、国内近、远期的需求量及主要销向预测。
㈡产品价格分析
1、国内矿产品价格现状。
2、矿产品价格稳定性及变化趋势。
三、矿产资源概况
㈠矿区总体概况
1、矿区总体规划情况。
2、矿区矿产资源概况。
3、该设计与矿区总体开发的关系。
㈡该设计项目的资源概况
1、矿床地质及构造特征。
2、矿床开采技术条件及水文地质条件。