A星算法matlab源码及详细注释
function astardemo
%ASTARDEMO Demonstration of ASTAR algorithm
%
% Copyright Bob L. Sturm, Ph. D., Assistant Professor
% Department of Architecture, Design and Media Technology
% formerly Medialogy
% Aalborg University i Ballerup
% formerly Aalborg University Copenhagen
% $Revision: 0.1 $ $Date: 2011 Jan. 15 18h24:24$
n = 20; % field size n x n tiles 20*20的界面
wallpercent = 0.45; % this percent of field is walls 45%的界面作为阻碍物（墙）% create the n x n FIELD with wallpercent walls containing movement costs,
% a starting position STARTPOSIND, a goal position GOALPOSIND, the costs
% A star will compute movement cost for each tile COSTCHART,
% and a matrix in which to store the pointers FIELDPOINTERS
[field, startposind, goalposind, costchart, fieldpointers] = ...
initializeField(n,wallpercent); %初始化界面
% initialize the OPEN and CLOSED sets and their costs
setOpen = [startposind]; setOpenCosts = [0]; setOpenHeuristics = [Inf];
setClosed = []; setClosedCosts = [];
movementdirections = {'R','L','D','U'};
% keep track of the number of iterations to exit gracefully if no solution counterIterations = 1;
% create figure so we can witness the magic
axishandle = createFigure(field,costchart,startposind,goalposind);
% as long as we have not found the goal or run out of spaces to explore
while ~max(ismember(setOpen,goalposind)) && ~isempty(setOpen) %ismember(A,B)返回与A同大小的矩阵，其中元素1表示A 中相应位置的元素在B中也出现，0则是没有出现
% for the element in OPEN with the smallest cost
[temp, ii] = min(setOpenCosts + setOpenHeuristics); %从OPEN表中选择花费最低的点temp,ii是其下标(也就是标号索引)
% find costs and heuristic of moving to neighbor spaces to goal
% in order 'R','L','D','U'
[costs,heuristics,posinds] = findFValue(setOpen(ii),setOpenCosts(ii), ...
field,goalposind,'euclidean'); %扩展temp的四个方向点，获得其坐标posinds，各个方向点的实际代价costs，启发代价heuristics
% put node in CLOSED and record its cost
setClosed = [setClosed; setOpen(ii)]; %将temp插入CLOSE表中setClosedCosts = [setClosedCosts; setOpenCosts(ii)]; %将temp的花费计入ClosedCosts
% update OPEN and their associated costs 更新OPEN表分为三种情况
if (ii > 1 && ii < length(setOpen)) %temp在OPEN表的中间，删除temp
setOpen = [setOpen(1:ii-1); setOpen(ii+1:end)];
setOpenCosts = [setOpenCosts(1:ii-1); setOpenCosts(ii+1:end)]; setOpenHeuristics = [setOpenHeuristics(1:ii-1); setOpenHeuristics(ii+1:end)];
elseif (ii == 1)
setOpen = setOpen(2:end); %temp是OPEN表的第一个元素，删除temp setOpenCosts = setOpenCosts(2:end);
setOpenHeuristics = setOpenHeuristics(2:end);
else %temp是OPEN表的最后一个元素，删除temp
setOpen = setOpen(1:end-1);
setOpenCosts = setOpenCosts(1:en
d-1);
setOpenHeuristics = setOpenHeuristics(1:end-1);
end
% for each of these neighbor spaces, assign costs and pointers;
% and if some are in the CLOSED set and their costs are smaller,
% update their costs and pointers
for jj=1:length(posinds) %对于扩展的四个方向的坐标
% if cost infinite, then it's a wall, so ignore
if ~isinf(costs(jj)) %如果此点的实际代价不为Inf,也就是没有遇到墙
% if node is not in OPEN or CLOSED then insert into costchart and
% movement pointers, and put node in OPEN
if ~max([setClosed; setOpen] == posinds(jj)) %如果此点不在OPEN表和CLOSE表中fieldpointers(posinds(jj)) = movementdirections(jj); %将此点的方向存在对应的fieldpointers中
costchart(posinds(jj)) = costs(jj); %将实际代价值存入对应的costchart中setOpen = [setOpen; posinds(jj)]; %将此点加入OPEN表中
setOpenCosts = [setOpenCosts; costs(jj)]; %更新OPEN表实际代价setOpenHeuristics = [setOpenHeuristics; heuristics(jj)]; %更新OPEN表启发代价% else node has already been seen, so check to see if we have
% found a better route to it.
elseif max(setOpen == posinds(jj)) %如果此点在OPEN表中
I = find(setOpen == posinds(jj)); %找到此点在OPEN表中的位置
% update if we have a better route
if setOpenCosts(I) > costs(jj) %如果在OPEN表中的此点实际代价比现在所得的大
costchart(setOpen(I)) = costs(jj); %将当前的代价存入costchart中，注意此点在costchart中的坐标与其自身坐标是一致的（setOpen(I)其实就是posinds(jj)），下同fieldpointers
setOpenCosts(I) = costs(jj); %更新OPEN表中的此点代价，注意此点在setOpenCosts中的坐标与在setOpen中是一致的，下同setOpenHeuristics setOpenHeuristics(I) = heuristics(jj); %更新OPEN表中的此点启发代价(窃以为这个是没有变的)
fieldpointers(setOpen(I)) = movementdirections(jj); %更新此点的方向
end
% else node has already been CLOSED, so check to see if we have
% found a better route to it.
else %如果此点在CLOSE表中，说明已经扩展过此点
% find relevant node in CLOSED
I = find(setClosed == posinds(jj));
% update if we have a better route
if setClosedCosts(I) > costs(jj) %如果在CLOSE表中的此点实际代价比现在所得的大（有一个问题，经过此点扩展的点还需要更新当前代价呢!!!）costchart(setClosed(I)) = costs(jj); %将当前的代价存入costchart中setClosedCosts(I) = costs(jj); %更新CLOSE表中的此点代价
fieldpointers(setClosed(I)) = movementdirections(jj); %更新此点的方向
end
end
end
end
if ise
mpty(setOpen) break; end %当OPEN表为空，代表可以经过的所有点已经查询完毕
set(axishandle,'CData',[costchart costchart(:,end); costchart(end,:) costchart(end,end)]);
% hack to make image look right
set(gca,'CLim',[0 1.1*max(costchart(find(costchart < Inf)))]); %CLim 将CData中的值与colormap对应起来: [cmin cmax] Color axis limits （不过不太明白为什么要*1.1）
drawnow; %cmin is the value of the data mapped to the first color in the colormap. cmax is the value of the data mapped to the last color in the colormap end
if max(ismember(setOpen,goalposind)) %当找到目标点时
disp('Solution found!'); %disp：Display array，disp(X)直接将矩阵显示出来，不显示其名字，如果X为string，就直接输出文字X
% now find the way back using FIELDPOINTERS, starting from goal position
p = findWayBack(goalposind,fieldpointers);
% plot final path
plot(p(:,2)+0.5,p(:,1)+0.5,'Color',0.2*ones(3,1),'LineWidth',4); drawnow;
elseif isempty(setOpen)
disp('No Solution!');
end
% end of the main function %%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p = findWayBack(goalposind,fieldpointers)
% This function will follow the pointers from the goal position to the
% starting position
n = length(fieldpointers); % length of the field
posind = goalposind;
% convert linear index into [row column]
[py,px] = ind2sub([n,n],posind);
% store initial position
p = [py px];
% until we are at the starting position
while ~strcmp(fieldpointers{posind},'S') %当查询到的点不是'S'起点时
switch fieldpointers{posind}
case 'L' % move left 如果获得该点的来源点方向为左时
px = px - 1;
case 'R' % move right
px = px + 1;
case 'U' % move up
py = py - 1;
case 'D' % move down
py = py + 1;
end
p = [p; py px];
% convert [row column] to linear index
posind = sub2ind([n n],py,px);
end
% end of this function %%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [cost,heuristic,posinds] = findFValue(posind,costsofar,field, ...
goalind,heuristicmethod)
% This function finds the movement COST for each tile surrounding POSIND in
% FIELD, returns their position indices POSINDS. They are ordered: right,
% left, down, up.
n = length(field); % length of the field
% convert linear index into [row column]
[currentpos(1) currentpos(2)] = ind2sub([n n],posind); %获得当前点的行列坐标，注意currentpos(1)是列坐标，currentpos(2)是行坐标
[goalpos(1) goalpos(2)] = ind2sub([n n],goalind); %获得目标点的行列坐标% places to store movement cost value and position
cost = Inf*ones(4,1); heuristic = Inf*ones(4,1); pos = ones(4,2);
% if we can look left, we move from the right 向左查询，那么就是从右边来newx = currentpos(2) - 1; n
ewy = currentpos(1);
if newx > 0 %如果没有到边界
pos(1,:) = [newy newx]; %获得新的坐标
switch lower(heuristicmethod)
case 'euclidean' %欧几里得距离(不像啊，亲)
heuristic(1) = abs(goalpos(2)-newx) + abs(goalpos(1)-newy); %heuristic(1)为启发函数计算的距离代价
case 'taxicab'
heuristic(1) = abs(goalpos(2)-newx) + abs(goalpos(1)-newy);
end
cost(1) = costsofar + field(newy,newx); %costsofar为之前花费的代价，field(newy,newx)为环境威胁代价，cost(1)为经过此方向点的真实代价
end
% if we can look right, we move from the left 向右查询
newx = currentpos(2) + 1; newy = currentpos(1);
if newx <= n
pos(2,:) = [newy newx];
switch lower(heuristicmethod)
case 'euclidean'
heuristic(2) = abs(goalpos(2)-newx) + abs(goalpos(1)-newy);
case 'taxicab'
heuristic(2) = abs(goalpos(2)-newx) + abs(goalpos(1)-newy);
end
cost(2) = costsofar + field(newy,newx);
end
% if we can look up, we move from down 向上查询
newx = currentpos(2); newy = currentpos(1)-1;
if newy > 0
pos(3,:) = [newy newx];
switch lower(heuristicmethod)
case 'euclidean'
heuristic(3) = abs(goalpos(2)-newx) + abs(goalpos(1)-newy);
case 'taxicab'
heuristic(3) = abs(goalpos(2)-newx) + abs(goalpos(1)-newy);
end
cost(3) = costsofar + field(newy,newx);
end
% if we can look down, we move from up 向下查询
newx = currentpos(2); newy = currentpos(1)+1;
if newy <= n
pos(4,:) = [newy newx];
switch lower(heuristicmethod)
case 'euclidean'
heuristic(4) = abs(goalpos(2)-newx) + abs(goalpos(1)-newy);
case 'taxicab'
heuristic(4) = abs(goalpos(2)-newx) + abs(goalpos(1)-newy);
end
cost(4) = costsofar + field(newy,newx);
end
% return [row column] to linear index
posinds = sub2ind([n n],pos(:,1),pos(:,2)); %posinds为此点扩展的四个方向上的坐标
% end of this function %%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化界面
function [field, startposind, goalposind, costchart, fieldpointers] = ... initializeField(n,wallpercent)
% This function will create a field with movement costs and walls, a start
% and goal position at random, a matrix in which the algorithm will store
% f values, and a cell matrix in which it will store pointers
% create the field and place walls with infinite cost 初始化界面和墙
field = ones(n,n) + 10*rand(n,n);
% field(ind2sub([n n],ceil(n^2.*rand(floor(n*n*wallpercent),1)))) = Inf; %floor(x)下取整，即舍去正小数至最近整数，ceil(x)上取整，即加入正小数至最近整数，Inf 代表正无穷
field(ceil(n^2.*rand(floor(n*n*wallpercent),1))) = Inf; %ind2sub是用来将线性坐标
(总体位置序号)
转为多维坐标(包含行列的坐标)的，发现其实不用转为多维坐标就可以，矩阵field可以访问线性坐标
% create random start position and goal position 随机选择行列作为起点与终点startposind = sub2ind([n,n],ceil(n.*rand),ceil(n.*rand)); %sub2ind用来将行列坐标转换为线性坐标，这里是必要的，因为如果把startposind设置成[x,y]的形式，访问field([x,y])的时候
goalposind = sub2ind([n,n],ceil(n.*rand),ceil(n.*rand)); %它并不是访问x行y列元素，而是访问线性坐标为x和y的两个元素
% force movement cost at start and goal positions to not be walls 将初始坐标设置为0，以免成为墙
field(startposind) = 0; field(goalposind) = 0;
% put not a numbers (NaN) in cost chart so A* knows where to look
costchart = NaN*ones(n,n); %costchart用来存储各个点的实际代价，NaN代表不是数据（不明确的操作）
% set the cost at the starting position to be 0
costchart(startposind) = 0; %起点的实际代价
% make fieldpointers as a cell array 生成n*n的元胞
fieldpointers = cell(n,n); %fieldpointers用来存储各个点的来源方向
% set the start pointer to be "S" for start, "G" for goal 起点设置为"S",终点设置为"G"
fieldpointers{startposind} = 'S'; fieldpointers{goalposind} = 'G'; % everywhere there is a wall, put a 0 so it is not considered 墙设置为0 fieldpointers(field == Inf) = {0}; %很好的方式，field == Inf 返回墙的位置，fieldpointers(field == Inf)设置相应的位置
% end of this function
%%%%%%%%%%%%%%%%%%%%
function axishandle = createFigure(field,costchart,startposind,goalposind)
% This function creates a pretty figure
% If there is no figure open, then create one
if isempty(gcbf) %gcbf是当前返回图像的句柄
f1 = figure('Position',[450 150 500 500],'Units','Normalized', ...
'MenuBar','none'); %这里的Position属性值为一个四元数组rect = [left, bottom, width, height]，第一、二个参数表示窗口位置，都是从屏幕的左下角计算的
%normalized —Units map the lower-left corner of the figure window to (0,0) and the upper-right corner to (1.0,1.0).
Caxes2 = axes('position', [0.01 0.01 0.98 0.98],'FontSize',12, ... 'FontName','Helvetica'); %position根据前面figure设置的
单位，in normalized units where (0,0) is the lower-left corner and (1.0,1.0) is the upper-right
else
% get the current figure, and clear it 获得当前图像并清空
gcf; cla;
end
n = length(field);
% plot field where walls are black, and everything else is white 0是黑色
field(field < Inf) = 0; %注意，虽然修改了field，但是这里的field属于局部变量，根本没有影响主函数中的field
pcolor([1:n+1],[1:n+1],[field field(:,end); field(end,:) f
ield(end,end)]); %多了一行一列
% set the colormap for the ploting the cost and looking really nice
cmap = flipud(colormap('jet')); %flipud用于反转矩阵colormap为产生jet类型的颜色表jet ranges from blue to red
% make first entry be white, and last be black
cmap(1,:) = zeros(3,1); cmap(end,:) = ones(3,1); %改变颜色表将尾色变为(0,0,0)是黑色，起色变为(1,1,1)是白色
% apply the colormap, but make red be closer to goal 红色是更接近目标的颜色colormap(flipud(cmap));
% keep the plot so we can plot over it
%********不用反转就可以*********%
%cmap = colormap('jet');
%cmap(1,:) = ones(3,1); cmap(end,:) = zeros(3,1);
%colormap(cmap);
%*******************************%
hold on;
% now plot the f values for all tiles evaluated
axishandle = pcolor([1:n+1],[1:n+1],[costchart costchart(:,end); costchart(end,:) costchart(end,end)]);
% plot goal as a yellow square, and start as a green circle
[goalposy,goalposx] = ind2sub([n,n],goalposind); %注意返回的列和行的位置[startposy,startposx] = ind2sub([n,n],startposind);
plot(goalposx+0.5,goalposy+0.5,'ys','MarkerSize',10,'Line Width',6); %加0.5是为了把坐标移到方块中央，'ys'中y表示yellow,s表示Square(方形)
plot(startposx+0.5,startposy+0.5,'go','MarkerSize',10,'Lin eWidth',6); %'go'中g表示green,o表示Circle(圆形)
% add a button so that can re-do the demonstration
uicontrol('Style','pushbutton','String','RE-DO& #39;, 'FontSize',12, ...
'Position', [1 1 60 40], 'Callback','astardemo');
% end of this function
function [R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=ACATSP(C,NC_max,m,Alpha,Beta ,Rho,Q) %%===================================================== ==================== %% ACATSP.m %% Ant Colony Algorithm for Traveling Salesman Problem %% ChengAihua,PLA Information Engineering University,ZhengZhou,China %% Email:aihuacheng@https://www.doczj.com/doc/666006711.html, %% All rights reserved %%------------------------------------------------------------------------- %% 主要符号说明 %% C n个城市的坐标，n×4的矩阵 %% NC_max 最大迭代次数 %% m 蚂蚁个数 %% Alpha 表征信息素重要程度的参数 %% Beta 表征启发式因子重要程度的参数 %% Rho 信息素蒸发系数 %% Q 信息素增加强度系数 %% R_best 各代最佳路线 %% L_best 各代最佳路线的长度 %%===================================================== ==================== %%第一步：变量初始化 n=size(C,1);%n表示问题的规模（城市个数） D=zeros(n,n);%D表示完全图的赋权邻接矩阵 for i=1:n for j=1:n if i~=j D(i,j)=max( ((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5,min(abs(C(i,3)-C(j,3)),144- abs(C(i,3)-C(j,3))) );%计算城市间距离 else D(i,j)=eps; end D(j,i)=D(i,j); end end Eta=1./D;%Eta为启发因子，这里设为距离的倒数 Tau=ones(n,n);%Tau为信息素矩阵 Tabu=zeros(m,n);%存储并记录路径的生成 NC=1;%迭代计数器 R_best=zeros(NC_max,n);%各代最佳路线
先新建一个主程序M文件ACATSP.m 代码如下: function [R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=ACATSP(C,NC_max,m,Alpha,Beta,Rho,Q) %%================================================== ======================= %% 主要符号说明 %% C n个城市的坐标，n×2的矩阵 %% NC_max 蚁群算法MATLAB程序最大迭代次数 %% m 蚂蚁个数 %% Alpha 表征信息素重要程度的参数 %% Beta 表征启发式因子重要程度的参数 %% Rho 信息素蒸发系数 %% Q 表示蚁群算法MATLAB程序信息素增加强度系数 %% R_best 各代最佳路线 %% L_best 各代最佳路线的长度 %%================================================== =======================
%% 蚁群算法MATLAB程序第一步：变量初始化 n=size(C,1);%n表示问题的规模（城市个数） D=zeros(n,n);%D表示完全图的赋权邻接矩阵 for i=1:n for j=1:n if i~=j D(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5; else D(i,j)=eps; % i = j 时不计算，应该为0，但后面的启发因子要取倒数，用eps（浮点相对精度）表示 end D(j,i)=D(i,j); %对称矩阵 end end Eta=1./D; %Eta为启发因子，这里设为距离的倒数 Tau=ones(n,n); %Tau为信息素矩阵 Tabu=zeros(m,n); %存储并记录路径的生成
第一章插值方法 1.1Lagrange插值 1.2逐步插值 1.3分段三次Hermite插值 1.4分段三次样条插值 第二章数值积分 2.1 Simpson公式 2.2 变步长梯形法 2.3 Romberg加速算法 2.4 三点Gauss公式 第三章常微分方程德差分方法 3.1 改进的Euler方法 3.2 四阶Runge-Kutta方法 3.3 二阶Adams预报校正系统 3.4 改进的四阶Adams预报校正系统 第四章方程求根 4.1 二分法 4.2 开方法 4.3 Newton下山法 4.4 快速弦截法 第五章线性方程组的迭代法 5.1 Jacobi迭代 5.2 Gauss-Seidel迭代 5.3 超松弛迭代 5.4 对称超松弛迭代 第六章线性方程组的直接法 6.1 追赶法 6.2 Cholesky方法 6.3 矩阵分解方法 6.4 Gauss列主元消去法
第一章插值方法 1.1Lagrange插值 计算Lagrange插值多项式在x=x0处的值. MATLAB文件：（文件名：Lagrange_eval.m）function [y0,N]= Lagrange_eval(X,Y,x0) %X,Y是已知插值点坐标 %x0是插值点 %y0是Lagrange插值多项式在x0处的值 %N是Lagrange插值函数的权系数 m=length(X); N=zeros(m,1); y0=0; for i=1:m N(i)=1; for j=1:m if j~=i; N(i)=N(i)*(x0-X(j))/(X(i)-X(j)); end end y0=y0+Y(i)*N(i); end 用法》X=[…];Y=[…]; 》x0= ; 》[y0,N]= Lagrange_eval(X,Y,x0) 1.2逐步插值 计算逐步插值多项式在x=x0处的值. MATLAB文件：（文件名：Neville_eval.m）function y0=Neville_eval(X,Y,x0) %X,Y是已知插值点坐标 %x0是插值点 %y0是Neville逐步插值多项式在x0处的值 m=length(X); P=zeros(m,1); P1=zeros(m,1); P=Y; for i=1:m P1=P; k=1; for j=i+1:m k=k+1;
蚁群算法的matlab源码，同时请指出为何不能优化到已知的最好解 % % % the procedure of ant colony algorithm for VRP % % % % % % % % % % % % %initialize the parameters of ant colony algorithms load data.txt; d=data(:,2:3); g=data(:,4); m=31; % 蚂蚁数 alpha=1; belta=4;% 决定tao和miu重要性的参数 lmda=0; rou=0.9; %衰减系数 q0=0.95; % 概率 tao0=1/(31*841.04);%初始信息素 Q=1;% 蚂蚁循环一周所释放的信息素 defined_phrm=15.0; % initial pheromone level value QV=100; % 车辆容量 vehicle_best=round(sum(g)/QV)+1; %所完成任务所需的最少车数V=40; % 计算两点的距离 for i=1:32; for j=1:32;
dist(i,j)=sqrt((d(i,1)-d(j,1))^2+(d(i,2)-d(j,2))^2); end; end; %给tao miu赋初值 for i=1:32; for j=1:32; if i~=j; %s(i,j)=dist(i,1)+dist(1,j)-dist(i,j); tao(i,j)=defined_phrm; miu(i,j)=1/dist(i,j); end; end; end; for k=1:32; for k=1:32; deltao(i,j)=0; end; end; best_cost=10000; for n_gen=1:50; print_head(n_gen); for i=1:m; %best_solution=[]; print_head2(i);
function [y,val]=QACStic load att48 att48; MAXIT=300; % 最大循环次数 NC=48; % 城市个数 tao=ones(48,48);% 初始时刻各边上的信息最为1 rho=0.2; % 挥发系数 alpha=1; beta=2; Q=100; mant=20; % 蚂蚁数量 iter=0; % 记录迭代次数 for i=1:NC % 计算各城市间的距离 for j=1:NC distance(i,j)=sqrt((att48(i,2)-att48(j,2))^2+(att48(i,3)-att48(j,3))^2); end end bestroute=zeros(1,48); % 用来记录最优路径 routelength=inf; % 用来记录当前找到的最优路径长度 % for i=1:mant % 确定各蚂蚁初始的位置 % end for ite=1:MAXIT for ka=1:mant %考查第K只蚂蚁 deltatao=zeros(48,48); % 第K只蚂蚁移动前各边上的信息增量为零 [routek,lengthk]=travel(distance,tao,alpha,beta); if lengthk 第8章 常微分方程初值问题数值解法 8.1 引言 8.2 欧拉方法 8.3 龙格-库塔方法 8.4 单步法的收敛性与稳定性 8.5 线性多步法 8.1 引 言 考虑一阶常微分方程的初值问题 00(,),[,],(). y f x y x a b y x y '=∈=（1.1） （1.2） 如果存在实数 ，使得 121212(,)(,).,R f x y f x y L y y y y -≤-?∈（1.3） 则称 关于 满足李普希茨(Lipschitz ）条件， 称为 的李普希茨常数（简称Lips.常数）. 0>L f y L f (参阅教材386页) 计算方法及MATLAB 实现 所谓数值解法，就是寻求解 在一系列离散节点 )(x y <<<<<+121n n x x x x 上的近似值 . ,,,,,121+n n y y y y 相邻两个节点的间距 称为步长. n n n x x h -=+1 如不特别说明，总是假定 为定数， ),2,1( ==i h h i 这时节点为 . ) ,2,1,0(0 =+=i nh x x n 初值问题(1.1),(1.2)的数值解法的基本特点是采取 “步进式”. 即求解过程顺着节点排列的次序一步一步地向前推进. 00(,),[,], (). y f x y x a b y x y '=∈= 描述这类算法，只要给出用已知信息 ,,,21--n n n y y y 计算 的递推公式. 1+n y 一类是计算 时只用到前一点的值 ，称为单步法. 1+n y n y 另一类是用到 前面 点的值 ， 1+n y k 11,,,+--k n n n y y y 称为 步法. k 其次，要研究公式的局部截断误差和阶，数值解 与 精确解 的误差估计及收敛性，还有递推公式的计算 稳定性等问题. n y )(n x y 首先对方程 离散化，建立求数值解的递推 公式. ),(y x f y =' 蚁群算法matlab精讲及仿真 4.1基本蚁群算法 4.1.1基本蚁群算法的原理 蚁群算法是上世纪90年代意大利学者M.Dorigo,v.Maneizz。等人提出来的,在越来越多的领域里得到广泛应用。蚁群算法，是一种模拟生物活动的智能算法，蚁群算法的运作机理来源于现实世界中蚂蚁的真实行为，该算法是由Marco Dorigo 首先提出并进行相关研究的，蚂蚁这种小生物，个体能力非常有限，但实际的活动中却可以搬动自己大几十倍的物体，其有序的合作能力可以与人类的集体完成浩大的工程非常相似，它们之前可以进行信息的交流，各自负责自己的任务，整个运作过程统一有序，在一只蚂蚁找食物的过程中，在自己走过的足迹上洒下某种物质，以传达信息给伙伴，吸引同伴向自己走过的路径上靠拢，当有一只蚂蚁找到食物后，它还可以沿着自己走过的路径返回，这样一来找到食物的蚂蚁走过的路径上信息传递物质的量就比较大，更多的蚂蚁就可能以更大的机率来选择这条路径，越来越多的蚂蚁都集中在这条路径上，蚂蚁就会成群结队在蚁窝与食物间的路径上工作。当然，信息传递物质会随着时间的推移而消失掉一部分，留下一部分，其含量是处于动态变化之中，起初，在没有蚂蚁找到食物的时候，其实所有从蚁窝出发的蚂蚁是保持一种随机的运动状态而进行食物搜索的，因此，这时，各蚂蚁间信息传递物质的参考其实是没有价值的，当有一只蚂蚁找到食物后，该蚂蚁一般就会向着出发地返回，这样，该蚂蚁来回一趟在自己的路径上留下的信息传递物质就相对较多，蚂蚁向着信息传递物质比较高的路径上运动，更多的蚂蚁就会选择找到食物的路径，而蚂蚁有时不一定向着信 息传递物质量高的路径走，可能搜索其它的路径。这样如果搜索到更短的路径后，蚂蚁又会往更短的路径上靠拢，最终多数蚂蚁在最短路径上工作。【基于蚁群算法和遗传算法的机器人路径规划研究】 该算法的特点： (1)自我组织能力，蚂蚁不需要知道整体环境信息，只需要得到自己周围的信息，并且通过信息传递物质来作用于周围的环境，根据其他蚂蚁的信息素来判断自己的路径。 (2)正反馈机制，蚂蚁在运动的过程中，收到其他蚂蚁的信息素影响，对于某路径上信息素越强的路径，其转向该路径的概率就越大，从而更容易使得蚁群寻找到最短的避障路径。 (3)易于与其他算法结合，现实中蚂蚁的工作过程简单，单位蚂蚁的任务也比较单一，因而蚁群算法的规则也比较简单，稳定性好，易于和其他算法结合使得避障路径规划效果更好。 (4)具有并行搜索能力探索过程彼此独立又相互影响，具备并行搜索能力，这样既可以保持解的多样性，又能够加速最优解的发现。 4.1.2 基本蚁群算法的生物仿真模型 a为蚂蚁所在洞穴,food为食物所在区,假设abde为一条路径,eadf为另外一条路径,蚂蚁走过后会留下信息素,5分钟后蚂蚁在两条路径上留下的信息素的量都为3,概率可以认为相同,而30分钟后baed路径上的信息素的量为60,明显大于eadf路径上的信息素的量。最终蚂蚁会完全选择abed这条最短路径,由此可见, 用MATLAB实现结构可靠度计算 口徐华…朝泽刚‘u刘勇‘21 。 (【l】中国地质大学(武汉工程学院湖北?武汉430074; 12】河海大学土木工程学院江苏?南京210098 摘要:Matlab提供了各种矩阵的运算和操作,其中包含结构可靠度计算中常用的各种数值计算方法工具箱,本文从基本原理和相关算例分析两方面,阐述利用Matlab,编制了计算结构可靠度Matlab程.序,使得Matlab-语言在可靠度计算中得到应用。 关键词:结构可靠度Matlab软件最优化法 中图分类号:TP39文献标识码:A文章编号:1007-3973(200902-095-Ol 1结构可靠度的计算方法 当川概率描述结构的可靠性时,计算结构可靠度就是计算结构在规定时问内、规定条件F结构能够完成预定功能的概率。 从简单到复杂或精确稃度的不同,先后提出的可靠度计算方法有一次二阶矩方法、二次二阶矩方法、蒙特卡洛方法以及其他方法。一次■阶矩方法又分为。I-心点法和验算点法,其中验算点法足H前可靠度分析最常川的方法。 2最优化方法计算可靠度指标数学模型 由结构111n个任意分布的独立随机变量一,x:…以表示的结构极限状态方程为:Z=g(■.托…t=0,采用R-F将非正念变量当罱正态化,得到等效正态分布的均值o:和标准差虹及可靠度指标B,由可靠度指标B的几何意义知。o;辟 开始时验算点未知,把6看成极限状态曲面上点P(■,爿:---37,的函数,通过优化求解,找到B最小值。求解可靠皮指标aJ以归结为以下约束优化模型: rain睁喜t华,2 s.,.Z=g(工i,x2’,…,工:=0 如极限状态方栉巾某个变最(X。可用其他变量表示,则上述模型jfIJ‘转化为无约束优化模型: 。。B!:手f生丛r+阻:坚:坠:盐尘}二剐 t∞oY?’【叫,J 3用MATLAB实现结构可靠度计算 3.1Matlab简介 Matlab是++种功能强、效率高、便.丁.进行科学和工程计算的交互式软件包,汇集了人量数学、统计、科学和工程所需的函数,MATI.AB具有编程简甲直观、用户界mf友善、开放性强等特点。将MATLAB用于蒙特卡罗法的一个显著优点是它拥有功能强大的随机数发生器指令。 3.2算例 3.2.I例:已知非线形极限状态方程z=g(t r'H=567f r-0.5H2=0’f、r服从正态分布。IIf=0.6,o r=0.0786;la|_ 2.18,o r_0.0654;H服从对数正态分布。u H= 3218,O。 =0.984。f、r、H相互独立,求可靠度指标B及验算点(,,r’,H‘。 解:先将H当量正念化:h=ln H服从正态分布,且 ,‘-““了:等专虿’=,。49?口二-、『五ir面_。。3 基于蚁群算法的机器人路径规划MATLAB源代码 基本思路是，使用离散化网格对带有障碍物的地图环境建模，将地图环境转化为邻接矩阵，最后使用蚁群算法寻找最短路径。 function [ROUTES,PL,Tau]=ACASPS(G,Tau,K,M,S,E,Alpha,Beta,Rho,Q) %% --------------------------------------------------------------- % ACASP.m % 基于蚁群算法的机器人路径规划 % GreenSim团队——专业级算法设计&代写程序 % 欢迎访问GreenSim团队主页→https://www.doczj.com/doc/666006711.html,/greensim %% --------------------------------------------------------------- % 输入参数列表 % G 地形图为01矩阵，如果为1表示障碍物 % Tau 初始信息素矩阵（认为前面的觅食活动中有残留的信息素） % K 迭代次数（指蚂蚁出动多少波） % M 蚂蚁个数（每一波蚂蚁有多少个） % S 起始点（最短路径的起始点） % E 终止点（最短路径的目的点） % Alpha 表征信息素重要程度的参数 % Beta 表征启发式因子重要程度的参数 % Rho 信息素蒸发系数 % Q 信息素增加强度系数 % % 输出参数列表 % ROUTES 每一代的每一只蚂蚁的爬行路线 % PL 每一代的每一只蚂蚁的爬行路线长度 % Tau 输出动态修正过的信息素 %% --------------------变量初始化---------------------------------- %load D=G2D(G); N=size(D,1);%N表示问题的规模（象素个数） MM=size(G,1); a=1;%小方格象素的边长 Ex=a*(mod(E,MM)-0.5);%终止点横坐标 if Ex==-0.5 Ex=MM-0.5; end Ey=a*(MM+0.5-ceil(E/MM));%终止点纵坐标 Eta=zeros(1,N);%启发式信息，取为至目标点的直线距离的倒数 %下面构造启发式信息矩阵 for i=1:N 作者：夏云木子 1、 >> syms re(x) re(y) re(z) >> input('计算相对误差：'),re(x)=10/1991,re(y)=0.0001/1.991,re(y)=0.0000001/0.0001991 所以可知re(y)最小，即y精度最高 2、 >> format short,A=sqrt(2) >> format short e,B=sqrt(2) >> format short g,C=sqrt(2) >> format long,D=sqrt(2) >> format long e,E=sqrt(2) >> format long g,F=sqrt(2) >> format bank,H=sqrt(2) >> format hex,I=sqrt(2) >> format +,J=sqrt(2) >> format,K=sqrt(2) 3、 >> syms A >> A=[sqrt(3) exp(7);sin(5) log(4)];vpa(pi*A,6) 4、1/6251-1/6252=1/6251*6252 5、（1）1/（1+3x）-(1-x)/(1+x)=x*(3*x-1)/[(1+3*x)*(1+x)] (2) sqrt(x+1/x)-sqrt(x-1/x)=2/x/[sqrt(x-1/x)+sqrt(x+1/x)] (3) log10(x1)-log(x2)=log10(x1/x2) (4) [1-cos(2*x)]/x =x^2/factorial(2)-x^4/factorial(4)+x^6/factorial(6)-… function [y,val]=QACS tic load att48 att48; MAXIT=300; % 最大循环次数 NC=48; % 城市个数 tao=ones(48,48);% 初始时刻各边上的信息最为1 rho=0.2; % 挥发系数 alpha=1; beta=2; Q=100; mant=20; % 蚂蚁数量 iter=0; % 记录迭代次数 for i=1:NC % 计算各城市间的距离 for j=1:NC distance(i,j)=sqrt((att48(i,2)-att48(j,2))^2+(att48(i,3)-att48(j,3))^2); end end bestroute=zeros(1,48); % 用来记录最优路径 routelength=inf; % 用来记录当前找到的最优路径长度 % for i=1:mant % 确定各蚂蚁初始的位置 % end for ite=1:MAXIT for ka=1:mant %考查第K只蚂蚁 deltatao=zeros(48,48); % 第K只蚂蚁移动前各边上的信息增量为零 [routek,lengthk]=travel(distance,tao,alpha,beta); if lengthk 1、Newdon迭代法求解非线性方程 function [x k t]=NewdonToEquation(f,df,x0,eps) %牛顿迭代法解线性方程 %[x k t]=NewdonToEquation(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:原函数，定义为内联函数 ?:函数的倒数，定义为内联函数 %x0:初始值 %eps:误差限 % %应用举例： %f=inline('x^3+4*x^2-10'); ?=inline('3*x^2+8*x'); %x=NewdonToEquation(f,df,1,0.5e-6) %[x k]=NewdonToEquation(f,df,1,0.5e-6) %[x k t]=NewdonToEquation(f,df,1,0.5e-6) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquation(f,df,1) if nargin==3 eps="0".5e-6; end tic; k=0; while 1 x="x0-f"(x0)./df(x0); k="k"+1; if abs(x-x0) < eps || k >30 break; end x0=x; end t=toc; if k >= 30 disp('迭代次数太多。'); x="0"; t="0"; end 2、Newdon迭代法求解非线性方程组 function y="NewdonF"(x) %牛顿迭代法解非线性方程组的测试函数 %定义是必须定义为列向量 y(1,1)=x(1).^2-10*x(1)+x(2).^2+8; y(2,1)=x(1).*x(2).^2+x(1)-10*x(2)+8; return; function y="NewdonDF"(x) %牛顿迭代法解非线性方程组的测试函数的导数 y(1,1)=2*x(1)-10; y(1,2)=2*x(2); y(2,1)=x(2).^+1; y(2,2)=2*x(1).*x(2)-10; return; 以上两个函数仅供下面程序的测试 function [x k t]=NewdonToEquations(f,df,x0,eps) %牛顿迭代法解非线性方程组 %[x k t]=NewdonToEquations(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:方程组（事先定义） ?:方程组的导数（事先定义） %x0:初始值 %eps:误差限 % %说明：由于虚参f和df的类型都是函数，使用前需要事先在当前目录下采用函数M文件定义% 另外在使用此函数求解非线性方程组时，需要在函数名前加符号“@”，如下所示 % %应用举例： %x0=[0,0];eps=0.5e-6; %x=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) 蚁群算法最短路径通用Matlab程序（附图） function [ROUTES,PL,Tau]=ACASP(G,Tau,K,M,S,E,Alpha,Beta,Rho,Q) %% --------------------------------------------------------------- % ACASP.m % 蚁群算法动态寻路算法 % ChengAihua,PLA Information Engineering University,ZhengZhou,China % Email:aihuacheng@https://www.doczj.com/doc/666006711.html, % All rights reserved %% --------------------------------------------------------------- % 输入参数列表 % G 地形图为01矩阵，如果为1表示障碍物 % Tau 初始信息素矩阵（认为前面的觅食活动中有残留的信息素） % K 迭代次数（指蚂蚁出动多少波） % M 蚂蚁个数（每一波蚂蚁有多少个） % S 起始点（最短路径的起始点） % E 终止点（最短路径的目的点） % Alpha 表征信息素重要程度的参数 % Beta 表征启发式因子重要程度的参数 % Rho 信息素蒸发系数 % Q 信息素增加强度系数 % % 输出参数列表 % ROUTES 每一代的每一只蚂蚁的爬行路线 % PL 每一代的每一只蚂蚁的爬行路线长度 % Tau 输出动态修正过的信息素 %% --------------------变量初始化---------------------------------- %load D=G2D(G); N=size(D,1);%N表示问题的规模（象素个数） MM=size(G,1); a=1;%小方格象素的边长 Ex=a*(mod(E,MM)-0.5);%终止点横坐标 if Ex==-0.5 Ex=MM-0.5; end Ey=a*(MM+0.5-ceil(E/MM));%终止点纵坐标 Eta=zeros(1,N);%启发式信息，取为至目标点的直线距离的倒数 %下面构造启发式信息矩阵 for i=1:N if ix==-0.5 信息系统分析与设计作业 层次分析法确定绩效评价权重在matlab中的实现 小组成员：孙高茹、王靖、李春梅、郭荣1 程序简要概述 编写程序一步实现评价指标特征值lam、特征向量w以及一致性比率CR的求解。 具体的操作步骤是：首先构造评价指标，用专家评定法对指标两两打分，构建比较矩阵，继而运用编写程序实现层次分析法在MATLAB中的应用。 通过编写MATLAB程序一步实现问题求解，可以简化权重计算方法与步骤，减少工作量，从而提高人力资源管理中绩效考核的科学化电算化。 2 程序在matlab中实现的具体步骤 function [w,lam,CR] = ccfx(A) %A为成对比较矩阵，返回值w为近似特征向量 % lam为近似最大特征值λmax，CR为一致性比率 n=length(A(:,1)); a=sum(A); B=A %用B代替A做计算 for j=1:n %将A的列向量归一化 B(:,j)=B(:,j)./a(j); end s=B(:,1); for j=2:n s=s+B(:,j); end c=sum(s);%计算近似最大特征值λmax w=s./c; d=A*w lam=1/n*sum((d./w)); CI=(lam-n)/(n-1);%一致性指标 RI=[0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45,1.49,1.51];%RI为随机一致 性指标 CR=CI/RI(n);%求一致性比率 if CR>0.1 disp('没有通过一致性检验'); else disp('通过一致性检验'); end end 3 案例应用 我们拟构建公司员工绩效评价分析权重，完整操作步骤如下： 3.1构建的评价指标体系 我们将影响员工绩效评定的指标因素分为：打卡、业绩、创新、态度与品德。 3.2专家打分，构建两两比较矩阵 A = 1.0000 0.5000 3.0000 4.0000 2.0000 1.0000 5.0000 3.0000 0.3333 0.2000 1.0000 2.0000 0.2500 0.3333 0.5000 1.0000 3.3在MATLAB中运用编写好的程序实现 直接在MATLAB命令窗口中输入 [w,lam,CR]=ccfx(A) 继而直接得出 d = 1.3035 2.0000 0.5145 0.3926 w = 0.3102 0.4691 0.1242 0.0966 lam =4.1687 计算方法作业（作者：夏云木子） 1、help linspace type linspace 2、a1=[5 12 47;13 41 2;9 6 71];a2=[12 9;6 15;7 21];B=a1*a2, C=a1(:,1:2).*a2, D=a1.^2, E=a1(:).^2 3、a1=[5 12 47;13 41 2;9 6 71];a2=[12 9;6 15;7 21];a1(4:5,1:3)=a2.';a1([4 5],:)=a1([5 4],:);b1=a1 c1=b1(4,1),c2=b1(5,3),D=b1(3:4,:)*a2 4、a1=[5 12 47;13 41 2;9 6 71]; E=eye(3,3); S = a1 + 5*a1' - E, S1=a1^3-rot90(a1)^2+6*E 5、a1=[5 12 47;13 41 2;9 6 71];s=5;A=s-a1,B=s*a1,C=s.*a1,D=s./a1,E=a1./s 6、c=[1 2 3 4;5 6 7 8;9 10 11 12;13 14 15 16];A=c^-4,B=(c^3)^-1,C=(3*c+5*c^-1)/5 7、a=[1 i 3;9i 2-i 8;7 4 8+i];A=a.' 8、abc=[-2.57 8.87;-0.57 3.2-5.5i];m1=sign(abc),m2=round(abc),m3=floor(abc) Sign为符号函数，round表示四舍五入取整，floor表示舍去小数部分取整 9、x=[1 4 3 2 0 8 10 5]';y=[8 0 0 4 2 1 9 11]';A=dot(x,y) 10、a=[3.82 5.71 9.62];b=[7.31 6.42 2.48];A=dot(a,b),B=cross(a,b) 11、P=[5 7 8 0 1];Pf=poly(P);Px=poly2str(Pf,'x') 12、P=[3 0 9 60 0 -90];K1=polyval(P,45),K2=polyval(P,-123),K3=polyval(P,579) 13、P1=[13 55 0 -17 9];P2=[63 0 26 -85 0 105];PP=conv(P1,P2);P1P2=poly2str(PP,'x'),[Q,r]=deconv(P2,P1) 蚁群算法最短路径通用Matlab程序 下面的程序是蚁群算法在最短路中的应用，稍加扩展即可应用于机器人路径规划 function [ROUTES,PL,Tau]=ACASP(G,Tau,K,M,S,E,Alpha,Beta,Rho,Q) %% ---------------------------------------------------------------% ACASP.m % 蚁群算法动态寻路算法 % ChengAihua,PLA Information Engineering University,ZhengZhou,China % Email:aihuacheng@https://www.doczj.com/doc/666006711.html, % All rights reserved %% ---------------------------------------------------------------% 输入参数列表 % G 地形图为01矩阵，如果为1表示障碍物 % Tau 初始信息素矩阵（认为前面的觅食活动中有残留的信息素） % K 迭代次数（指蚂蚁出动多少波） % M 蚂蚁个数（每一波蚂蚁有多少个） % S 起始点（最短路径的起始点） % E 终止点（最短路径的目的点） % Alpha 表征信息素重要程度的参数 % Beta 表征启发式因子重要程度的参数 % Rho 信息素蒸发系数 % Q 信息素增加强度系数 % % 输出参数列表 % ROUTES 每一代的每一只蚂蚁的爬行路线 % PL 每一代的每一只蚂蚁的爬行路线长度 % Tau 输出动态修正过的信息素 %% --------------------变量初始化---------------------------------- %load D=G2D(G); N=size(D,1);%N表示问题的规模（象素个数） MM=size(G,1); a=1;%小方格象素的边长 Ex=a*(mod(E,MM)-0.5);%终止点横坐标 if Ex==-0.5 Ex=MM-0.5; end Ey=a*(MM+0.5-ceil(E/MM));%终止点纵坐标 Eta=zeros(1,N);%启发式信息，取为至目标点的直线距离的倒数 %下面构造启发式信息矩阵 for i=1:N if ix==-0.5 ix=MM-0.5; 基于蚁群算法的PID控制参数优化Matlab源码 (2009-07-26 12:31:02) 除了蚁群算法，可用于PID参数优化的智能算法还有很多，比如遗传算法、模拟退火算法、粒子群算法、人工鱼群算法，等等。 function [BESTX,BESTY,ALLX,ALLY]=ACOUCP (K,N,Rho,Q,Lambda,LB,UB,Num,Den,Delay,ts,StepNum,SigType,PIDLB,PIDUB) %% 此函数实现蚁群算法，用于PID控制参数优化 % GreenSim团队原创作品，转载请注明 % Email:greensim@https://www.doczj.com/doc/666006711.html, % GreenSim团队主页：https://www.doczj.com/doc/666006711.html,/greensim % [color=red]欢迎访问GreenSim——算法仿真团队→[url=https://www.doczj.com/doc/666006711.html,/greensim] https://www.doczj.com/doc/666006711.html,/greensim[/url][/color] %% 输入参数列表 % K 迭代次数 % N 蚁群规模 % Rho 信息素蒸发系数，取值0～1之间，推荐取值0.7～0.95 % Q 信息素增加强度，大于0，推荐取值1左右 % Lambda 蚂蚁爬行速度，取值0～1之间，推荐取值0.1～0.5 % LB 决策变量的下界，M×1的向量 % UB 决策变量的上界，M×1的向量 % Num 被控制对象传递函数的分子系数向量 % Den 被控制对象传递函数的分母系数向量 % Delay 时间延迟 % ts 仿真时间步长 % StepNum 仿真总步数 % SigType 信号类型，1为阶跃信号，2为方波信号，3为正弦波信号 % PIDLB PID控制输出信号限幅的下限 % PIDUB PID控制输出信号限幅的上限 %% 输出参数列表 % BESTX K×1细胞结构，每一个元素是M×1向量，记录每一代的最优蚂蚁 % BESTY K×1矩阵，记录每一代的最优蚂蚁的评价函数值 % ALLX K×1细胞结构，每一个元素是M×N矩阵，记录每一代蚂蚁的位置 % ALLY K×N矩阵，记录每一代蚂蚁的评价函数值 function [x0,k]=bisect1(fun1,a,b,ep) if nargin<4 ep=1e-5; end fa=feval(fun1,a); fb=feval(fun1,b); if fa*fb>0 x0=[fa,fb]; k=0; return; end k=1; while abs(b-a)/2>ep x=(a+b)/2; fx=feval(fun1,x); if fx*fa<0 b=x; fb=fx; else a=x; fa=fx; k=k+1; end end x0=(a+b)/2; >> fun1=inline('x^3-x-1'); >> [x0,k]=bisect1(fun1,1.3,1.4,1e-4) x0 = 1.3247 k = 7 >> 简单迭代法 function [x0,k]=iterate1(fun1,x0,ep,N) if nargin<4 N=500; end if nargin<3 ep=1e-5; end x=x0; x0=x+2*ep; while abs(x-x0)>ep & k0计算方法及MATLAB实现简明讲义课件PPS8-1欧拉龙格法
matlab蚁群算法精讲及仿真图
用MATLAB实现结构可靠度计算.
基于蚁群算法的MATLAB实现
计算方法及其MATLAB实现第二章作业
(完整版)蚁群算法matlab程序实例整理
matlab用于计算方法的源程序
蚁群算法最短路径通用Matlab程序(附图)
层次分析法计算权重在matlab中的实现
计算方法及其MATLAB实现第一章作业
蚁群算法最短路径matlab程序
基于蚁群算法的PID控制参数优化Matlab源码
数值计算方法matlab程序
相关主题
文本预览