2.虚数矩阵
C_real=real(C)
C_imag=imag(C)
C_magnitude=abs(C)
3.画出震荡曲线y=exp(-t/3)*sin(3*t)
t=0:pi/50:4*pi;
y0=exp(-t/3);
y=exp(-t/3) .*sin(3*t);
plot(t,y,'-r',t,y0,':b',t,-y0,':b')
grid
4.画出三锥曲面
4.0
clear;
x=-8:0.5:8; %定义自变量x的一位刻度向量
y=x'; %定义自变量y的一位刻度向量
X=ones(size(y))*x; %计算自变量平面上取值点x坐标的二维数组
Y=y*ones(size(x)); %计算自变量平面上取值点y坐标的二维数组
R=sqrt(X.^2+Y.^2)+eps; %计算中间变量R=x,y平方和的平方根
Z=sin(R)./R; %计算与自变量二维数组相应的函数值z=sinR/R
mesh(Z); %绘制三维网格图
colormap(hot)
4.1 Schaffer函数
N=10;
n=0.1;
M=2*N/n+1;
x=-N:n:N; %定义自变量x的一位刻度向量
y=-N:n:N; %定义自变量x的一位刻度向量
[X,Y]=meshgrid(x,y);
R=sqrt(X.^2+Y.^2)+eps; %计算中间变量R=x,y平方和的平方根
A=sin(R).^2;
B=(1+0.001*R.^2)^2;
z=0.5+(A-0.5)./B;
mesh(z); %绘制三维网格图
colormap(hot)
4.2 Shubert函数
N=10;
n=0.1;
M=2*N/n+1;
x=-N:n:N; %定义自变量x的一位刻度向量
y=-N:n:N; %定义自变量x的一位刻度向量
[X,Y]=meshgrid(x,y);
A=zeros(M,M);B=zeros(M,M);
for i=1:5
R=(i+1)*X+i;
S=(i+1)*Y+i;
A=A+i*cos(R);
B=B+i*cos(S);
end
z=A.*B;
mesh(z); %绘制三维网格图
colormap(hot)
4.3 Griewank函数
N=10;
n=0.1;
M=2*N/n+1;
x=-N:n:N; %定义自变量x的一位刻度向量
y=-N:n:N; %定义自变量x的一位刻度向量
[X,Y]=meshgrid(x,y);
z=X.^2/4000+Y.^2/4000+cos(X).*cos(Y/2)+1;
mesh(z); %绘制三维网格图
colormap(hot)
4.4 Rastrigrin函数
N=2.12;
n=0.05;
M=2*N/n+1;
x=-N:n:N; %定义自变量x的一位刻度向量
y=-N:n:N; %定义自变量x的一位刻度向量
[X,Y]=meshgrid(x,y);
z=(X.^2-10*cos(2*pi*X)+10)+(Y.^2-10*cos(2*pi*Y)+10);
mesh(z);
colormap(hot)
5.数据存取
mkdir('c:\\','my_dir'); %在C盘上创建目录my_dir
cd c:\\my_dir %使c:\\my_dir成为当前目录
save saf X Y Z %选择内存中的X、Y、Z保存到saf.mat
dir %显示目录上的文件
%清空内存，装载变量Z
clear %清空内存
load saf Z %把变量Z读入内存
who %检查内存中的变量
第二节
1.矩阵的输入：
Time=[1 2 3 4 5 6 7 8 9 10 11 12 13]
X_Da=[2.32 3.43;4.37 5.98]
vect_a=[1 2 3 4 5]
Matrix_B=[1 2 3;
2 3 4;3 4 5]
Null_M=[] %生成空矩阵
2、复数矩阵
a=2.7;b=13/25;
C=[1,2*a+i*b,b*sqrt(a);sin(pi/4),a+5*b,3.5+1]
或者
R=[1 2 3;4 5 6],M=[11 12 13;14 15 16]
CN=R+i*M
3、sym函数——定义符号矩阵
sym(符号量1；符号量2；符号量n）
sym_matrix=sym('[A B C;NIHAO,Welcome,
to beijing],')
sym_digits=sym('[1 2 3;a b c;sin(x) cos(y) tan(z)]')
clc;
clear;
sym_matrix=sym('[a,b,c;Jack,Help,Me!;NO,WAY,!]')
把数值矩阵定义为符号矩阵
a=[1/3 sqrt(2) 3.4234;exp(0.23) log(29) 23^(-11.23)]
s=sym(a)
4、cat函数
A=cat(n,A1,A2,···,Am) %创建数组
A1=[1,2,3;4,5,6;7,8,9];A2=A1';A3=A1-A2;
A4=cat(3,A1,A2,A3)
5、zeros()
zeros(n)\\zeros(m,n)\\zeros([m,n])\\zeros(size(A))
n=6;
A=zeros(n)
6、eye函数，单位矩阵生成
Y=eye(n)、eye（m,n)、eye(size(A))
n=4;
m=5;
Y1=eye(n)
Y2=eye(m,n)
Y3=eye(size(Y1))
7、ones函数，生成全1阵
Y=ones(n)、ones(m,n)、ones(size(A))等
n=4;
m=6;
X1=ones(m,n)
d1=2;d2=3;d3=4;
X2=ones(d1,d2,d3)
8、rand函数，生成均匀分布随机矩阵
Y=rand(n)\
and(m,n)\
and(size(A))等
s=rand('state') 产生包括均匀发生器当前状态的35个元素的向量。
rand('state',s) 使状态重置为S.
rand('state',0) 重置发生器到初始状态。
rand('state',j) 对整数j重置发生器到第j个状态。
rand('state',sum(100*clock)) 每次重置到不同的状态。
R=rand(3,4)
a=10;b=20;
x=a+(b-a)*rand(4)
9、randn函数，生成正态随机矩阵
Y=randn(n) %生成n*n正态分布随机矩阵
等等同上
randn %无变量输入时只产生一个正态分布随机数
s=randn('state') 产生包括均匀发生器当前状态的2个元素的向量。
rand('state',s) 使状态重置为S.
rand('state',0) 重置发生器到初始状态。
rand('state',j) 对整数j重置发生器到第j个状态。
rand('state',sum(100*clock)) 每次重置到不同的状态。
产生均值为0.6，方差为0.1的4阶矩阵
mu=0.6;sigma=0.1;
x=mu+sqrt(sigma)*randn(4)
10、randperm函数，产生随机序列
p=randperm(n) %产生1~n之间整数的随机排列
randperm(6)
11、linspace函数，线性等分向量的生成
y=linspace(a,b） %在（a,b)上产生100个线性等分点
y=linspace(a,b,n) %在（a,b)上产生 n 个线性等分点
在1-200之间生成100个\\20个等分点
a=1;b=200;
y=linspace(a,b)
z=linspace(a,b,100)
a=0;b=200;
y=linspace(a,b)
z=linspace(a,b,100)
12、logspace函数，对数等分向量的生成
y=logspace(a,b） %在（10^a,10^b)上产生50个对数等分点
y=l0gspace(a,b,n) %在（10^a,10^b)上产生n 个对数等分点
在10^2-10^4之间生成50个\\30个等分点
a=2;b=4;
x1=logspace(a,b)
x2=logspace(a,b,30)
13、compan函数，生成友矩阵（求函数的根）
A=compan(u)
求（x-1)（x+2)(x-3)、x^3-8*x+13的友矩阵和根
u=[1,0,-8,13];
A=compan(u)
eig(A)
u=[1,-2,-5,6];
A=compan(u)
eig(A)
第一次实现
简单实现
han1=[83.0 79.8 78.1 85.1 86.6 88.2 90.3 86.7 93.3 92.5 90.9 96.9
101.7 85.1 87.8 91.6 93.4 94.5 97.4 99.5 104.2 102.3 101.0 123.5
92.2 114.0 93.3 101.0 103.5 105.2 109.5 109.2 109.6 111.2 121.7 131.3
105.0 125.7 106.6 116.0 117.6 118.0 121.7 118.7 120.2 127.8 121.8 121.9
139.3 129.5 122.5 124.5 135.7 130.8 138.7 133.7 136.8 138.9 129.6 133.7
137.5 135.3 133.0 133.4 142.8 141.6 142.9 147.3 159.6 162.1 153.5 155.9
163.2 159.7 158.4 145.2 124.0 144.1 157.0 162.6 171.8 180.7 173.5 176.5];
han1(end,:)=[];m=size(han1,2); %令最后一行为空、求数列的列数
x0=mean(han1,2); %求han1的每行的平均数，即每年的均值
x1=cumsum(x0); %x1(k)=x0(k-1)+x0(k)
alpha=0.4;n=length(x0); %求数列z1
z1=alpha*x1(2:n)+(1-alpha)*x1(1:n-1) %求数列z1
Y=x0(2:n);B=[-z1,ones(n-1,1)];
ab=B\\Y %求出参数a、b
k=6; %此时求的是第k+1项预测值
x7hat=(x0(1)-ab(2)/ab(1))*(exp(-ab(1)*k)-exp(-ab(1)*(k-1)))%直接求第七年的预测值（这里是平均值）
z=m*x7hat %求第七年的总值
u=sum(han1)/sum(sum(han1)) %求每月比例=（各年同月之和）/（各年总和）
v=z*u %求个月预测值
例2：较复杂
renkou
renkou1=rk(:,1); %renkou数列第一列，即年末常住人口数
renkou2=rk(:,2); %renkou数列第二列，即户籍人口
renkou3=rk(:,3); %renkou数列第三列，即非户籍人口
x0=renkou2'; %赋予初始值，拷贝，并转置
n=length(x0); %下标
lamda=x0(1:n-1)./x0(2:n) %求数列的级比
range=minmax(lamda) %级比的范围
x1=cumsum(x0) %求累积和x1
for i=2:n
z(i)=0.5*(x1(i)+x1(i-1));
end %当α=0.5时的z1
B=[-z(2:n)',ones(n-1,1)]; %令数列
Y=x0(2:n)'; %令数列
u=B\\Y %求参数向量u=(a,b)
x=dsolve('Dx+a*x=b','x(0)=x0'); %求解常微分，并令x(1)=x0
x=subs(x,{'a','b','x0'},{u(1),u(2),x1(1)}); %进行符号替换
yuce1=subs(x,'t',[0:n-1]); %把自变量由连续的改为整数
digits(6),y=vpa(x) %将x保留5位小数
yuce=[x0(1),diff(yuce1)]
epsilon=x0-yuce %计算残差
delta=abs(epsilon./x0) %计算相对误差
rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda %计算级比偏差值
clc,clear
han1=[83.0 79.8 78.1 85.1 86.6 88.2 90.3 86.7 93.3 92.5 90.9 96.9
101.7 85.1 87.8 91.6 93.4 94.5 97.4 99.5 104.2 102.3 101.0 123.5
92.2 114.0 93.3 101.0 103.5 105.2 1
09.5 109.2 109.6 111.2 121.7 131.3
105.0 125.7 106.6 116.0 117.6 118.0 121.7 118.7 120.2 127.8 121.8 121.9
139.3 129.5 122.5 124.5 135.7 130.8 138.7 133.7 136.8 138.9 129.6 133.7
137.5 135.3 133.0 133.4 142.8 141.6 142.9 147.3 159.6 162.1 153.5 155.9
163.2 159.7 158.4 145.2 124.0 144.1 157.0 162.6 171.8 180.7 173.5 176.5];
han1(end,:)=[];m=size(han1,2); %令最后一行为空、求数列的列数
x0=mean(han1,2); %求han1的每行的平均数，即每年的均值
x1=cumsum(x0); %x1(k)=x0(k-1)+x0(k)
alpha=0.4;n=length(x0); %求数列z1
z1=alpha*x1(2:n)+(1-alpha)*x1(1:n-1) %求数列z1
Y=x0(2:n);B=[-z1,ones(n-1,1)];
ab=B\\Y %求出参数a、b
k=6;
x7hat=(x0(1)-ab(2)/ab(1))*(exp(-ab(1)*k)-exp(-ab(1)*(k-1)))%直接求第七年的预测值（这里是平均值）
z=m*x7hat %求第七年的总值
u=sum(han1)/sum(sum(han1)) %求每月比例=（各年同月之和）/（各年总和）
v=z*u %求个月预测值
renkou
renkou1=rk(:,1); %renkou数列第一列，即年末常住人口数
renkou2=rk(:,2); %renkou数列第二列，即户籍人口
renkou3=rk(:,3); %renkou数列第三列，即非户籍人口
x0=renkou2'; %赋予初始值，拷贝，并转置
n=length(x0); %下标
lamda=x0(1:n-1)./x0(2:n) %求数列的级比
range=minmax(lamda) %级比的范围
x1=cumsum(x0) %求累积和x1
for i=2:n
z(i)=0.5*(x1(i)+x1(i-1));
end %当α=0.5时的z1
B=[-z(2:n)',ones(n-1,1)]; %令数列
Y=x0(2:n)'; %令数列
u=B\\Y %求参数向量u=(a,b)
x=dsolve('Dx+a*x=b','x(0)=x0'); %求解常微分，并令x(1)=x0
x=subs(x,{'a','b','x0'},{u(1),u(2),x1(1)}); %进行符号替换
yuce1=subs(x,'t',[0:n-1]); %把自变量由连续的改为整数
digits(6),y=vpa(x) %将x保留5位小数
yuce(1)=x0(1);
for i=2:n
yuce(i)=yuce1(i)-yuce1(i-1);
end
yuce=vpa(yuce)
epsilon=x0-yuce %计算残差
delta=abs(epsilon./x0) %计算相对误差
rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda %计算级比偏差值
2、遗传算法
第一次尝试：简单版
fplot('variable.*variable',[0,31]); %画出函数varia
ble.*variable在?[0,31]上的图像
NIND=4; %种群个数目
MAXGEN=10; %最大遗代数
PRECI=5; %变量的二进制数
GGAP=1; %代沟
trace=zeros(2,MAXGEN); %寻优结果的初始值
FieldD=[5;0;31;1;0;1;1]; %区域描述器
Chrom=crtbp(NIND,PRECI); %初始种群
gen=0; %计数器
variable=bs2rv(Chrom,FieldD); %计算初始种群的十进制转换
ObjV=variable.*variable; %计算目标函数值
while gen SelCh=select('sus',Chrom,FitnV,GGAP); %选择 SelCh=recombin('xovsp',SelCh,1); %重组，100% SelCh=mut(SelCh); %变异 variable=bs2rv(SelCh,FieldD); %子代个体的十进制换 ObjVSel=variable.*variable; %计算子代的目标函数值 [Chrom ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjVSel); %重插入子代的新种群 variable=bs2rv(Chrom,FieldD); %子代个体的十进制转换 gen=gen+1; %代计数器加1 %输出最优解及其序号，并在目标函数图像中标出，Y为最优解,I为种群的序号 [Y,I]=max(ObjV); hold on; plot(variable(I),Y,'bo'); trace(1,gen)=max(ObjV); %遗传算法性能跟踪 trace(2,gen)=sum(ObjV)/length(ObjV); end variable=bs2rv(Chrom,FieldD); %最优个体的十进制转换 grid,hold on; plot(variable,ObjV,'b*'); figure(2); plot(trace(1,:)); hold on; plot(trace(2,:),'-.');grid legend('解的变化','种群均值的变化') 第二次尝试，复杂类 优化函数使用简单的遗传算法与最优保留（一次更新,一次交叉，一次变异） Max f(x1,x2)=100*(x1*x1-x2).^2+(1-x1).^2; -2.0480<=x1,x2<=2.0480 format long; %设定数据显示格式 %初始化参数 T=100; %仿真代数 N=80; %群体规模 pm=0.05;pc=0.8; %交叉变异概率 umax=2.048;umin=-2.048; %参数取值范围 L=10; %单个参数字串长度，总编码长度2L bval=round(rand(N,2*L)); %初始种群 bestv=-inf; %最优适应度初值，即无穷大 %迭代开始 for ii=1:T %解码，计算适应度 for i=1:N %求第ii代时所有样本的目标函数值obj、及对应x1,x1即xx( ， , ) y1=0;y2=0; for j=1:1:L y1=y1+bval(i,L-j+1)*2^(j-1); %把第i行前10个元素的二进制化为十进制y1 end x1=(umax-umin)*y1/(2^L-1)+umin; %计算y1对应的x1，份数比例*区间大小(最大减去最小)+最小值(负值) for j=1:1:L y2=y2+bval(i,2*L-j+1)*2^(j-1); %把第i行后10个元素的二进制化为十进制y2 end x2=(umax-umin)*y2/(2^L-1)+umin; %计算y2对应的x2，份数比例*区间大小(最大减去最小)+最小值(负值) obj(i)=100*(x1*x1-x2).^2+(1-x1).^2; %目标函数，其实也是指标函数值 xx(i,:)=[x1,x2]; %第i行xx记录x1,x1 end func=obj; %目标函数转换为适应度函数 p=func./sum(func); %每个样本适应度值占总数的百分比 q=cumsum(p); %比例累加值 [fmax,indmax]=max(func); %求当代最佳个体，fmax显示最大值，indmax表示最大值对应位置，即i if fmax>=bestv %每一代计算一次，留下目前为止的最大适应度值 bestv=fmax; %到目前为止最优适应度值 bvalxx=bval(indmax,:); %到目前为止最佳位串（初始种群的二进制串） optxx=xx(indmax,:); %到目前为止最优参数（即最优个体对应的x1,x2） end Bfit1(ii)=bestv; % 存储每代的最优适应度 %%%%遗传操作开始 %轮盘赌选择 for i=1:(N-1) %随机产生N-1个数字，好处对应的样本二进制串 r=rand; %产生随机数r tmp=find(r<=q); newbval(i,:)=bval(tmp(1),:); %找到r对应的样本 end newbval(N,:)=bvalxx; %最优保留 bval=newbval; %更新样本，完成选择--复制 %单点交叉 for i=1:2:(N-1) %每两个样本进行一次交叉 cc=rand; %随机生成cc if cc ch=bval(i,:); bval(i,point+1:2*L)=bval(i+1,point+1:2*L); bval(i+1,point+1:2*L)=ch(1,point+1:2*L);%实现样本二进制串i和i+1后point+1到2L的交叉 end end bval(N,:)=bvalxx; %最优保留,当上面的i=79时，会把最后一个改变，现在再恢复最大 %位 点变异 mm=rand(N,2*L) bval(mm)=1-bval(mm); %与稀疏矩阵mm非零位对应位置突变 end %输出 plot(Bfit1); %绘制最优适应度进化曲线 bestv %输出最优适应度值 optxx %输出最优参数 第三个 sjs %加载敌方 100 个目标的数据 d1=[70,40]; sj0=[d1;sj;d1]; %距离矩阵 d sj=sj0*pi/180; d=zeros(102); for i=1:101 for j=i+1:102 temp=cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2)); d(i,j)=6370*acos(temp); end end d=d+d';L=102;w=50;dai=100; %通过改良圈算法选取优良父代 A for k=1:w c=randperm(100); c1=[1,c+1,102]; flag=1; while flag>0 flag=0; for m=1:L-3 for n=m+2:L-1 if d(c1(m),c1(n))+d(c1(m+1),c1(n+1)) c1(m+1:n)=c1(n:-1:m+1); end end end end J(k,c1)=1:102; end J=J/102; J(:,1)=0;J(:,102)=1; rand('state',sum(clock)); %遗传算法实现过程 A=J; for k=1:dai %产生 0～1间随机数列进行编码 B=A; c=randperm(w); %交配产生子代 B for i=1:2:w F=2+floor(100*rand(1)); temp=B(c(i),F:102); B(c(i),F:102)=B(c(i+1),F:102); B(c(i+1),F:102)=temp; end %变异产生子代 C by=find(rand(1,w)<0.1); if length(by)==0 by=floor(w*rand(1))+1; end C=A(by,:); L3=length(by); for j=1:L3 bw=2+floor(100*rand(1,3)); bw=sort(bw); C(j,:)=C(j,[1:bw(1)-1,bw(2)+1:bw(3),bw(1):bw(2),bw(3)+1:102]); end G=[A;B;C]; TL=size(G,1); %在父代和子代中选择优良品种作为新的父代 [dd,IX]=sort(G,2);temp(1:TL)=0; for j=1:TL for i=1:101 temp(j)=temp(j)+d(IX(j,i),IX(j,i+1)); end end [DZ,IZ]=sort(temp); A=G(IZ(1:w),:); end path=IX(IZ(1),:) long=DZ(1) xx=sj0(path,1);yy=sj0(path,2); plot(xx,yy,'-o') 3、粒子群算法 function DrawRastrigin() %绘制Rastrigin函数图形 x=[-5:0.05:5]; y=x; [X,Y]=meshgrid(x,y); [row,col]=size(X); for l=1:col for h=1:row z(h,l)=Rastrigin([X(h,l),Y(h,l)]); end end surf(X,Y,z); shading interp 第二次尝试：复杂类 %声明参数的优化 max_iterations = 1000; no_of_particles = 50; dimensions = 1; delta_min = -0.003; delta_max = 0.003; c1 = 1.3; c2 = 1.3; %初始化粒子及其速度分量 for count_x = 1:no_of_particles for count_y = 1:dimensions particle_position(count_x,count_y) = rand*10; particle_velocity(count_x,count_y) = rand; p_best(count_x,count_y) = particle_position(count_x, x,count_y); end end %initialize the p_best_fitness array for count = 1:no_of_particles p_best_fitness(count) = -1000; end %particle_position %particle_velocity %main particle swrm routine for count = 1:max_iterations %find the fitness of each particle %change fitness function as per equation requiresd and dimensions for count_x = 1:no_of_particles %x = particle_position(count_x,1); %y = particle_position(count_x,2); %z = particle_position(count_x,3); %soln = x^2 - 3*y*x + z; %x = particle_position(count_x); %soln = x^2-2*x+1; x = particle_position(count_x); soln = x-7; if soln~=0 current_fitness(count_x) = 1/abs(soln); else current_fitness =1000; end end %decide on p_best etc for each particle for count_x = 1:no_of_particles if current_fitness(count_x) >p_best_fitness(count_x) p_best_fitness(count_x) = current_fitness(count_x); for count_y = 1:dimensions p_best(count_x,count_y) = particle_position(count_x,count_y); end end end %decide on the global best among all the particles [g_best_val,g_best_index] = max(current_fitness); %g_best contains the position of teh global best for count_y = 1:dimensions g_best(count_y) = particle_position(g_best_index,count_y); end %update the position and velocity compponents for count_x = 1:no_of_particles for count_y = 1:dimensions p_current(count_y) = particle_position(count_x,count_y); end for count_y = 1:dimensions particle_velocity(count_y) = particle_velocity(count_y) + c1*rand*(p_best(count_y)-p_current(count_y)) + c2*rand*(g_best(count_y)-p_current(count_y)); particle_positon(count_x,count_y) = p_current(count_y) +particle_velocity(count_y); end end end g_best current_fitness(g_best_index) clear all, clc % pso example iter = 1000; % number of algorithm iterations np = 2; % number of model parameters ns = 10; % number of sets of model parameters Wmax = 0.9; % maximum inertial weight Wmin = 0.4; % minimum inertial weight c1 = 2.0; % parameter in PSO methodology c2 = 2.0; % parameter in PSO methodology Pmax = [10 10]; % maximum model parameter value Pmin = [-10 -10]; % minimum model parameter value Vmax = [1 1]; % maximum change in model parameter Vmin = [-1 -1]; % minimum change in model parameter modelparameters(1:np,1:ns) = 0; % set all model parameter estimates for all model parameter sets to zero modelparameterchanges(1:np,1:ns) = 0; % set all change in model parameter estimates for all model parameter sets to zero bestmodelparameters(1:np,1:ns) = 0; % set best model parameter estimates for all model parameter sets to zero setbestcostfunction(1:ns) = 1e6; % set best cost function of each model parameter set to a large number globalbestparameters(1:np) = 0; % set best model parameter values for all model parameter sets to zero bestparameters = globalbestparameters'; % best model parameter values for all model parameter sets (to plot) globalbestcostfunction = 1e6; % set best cost function for all model parameter sets to a large number i = 0; % indicates ith algorithm iteration j = 0; % indicates jth set of model parameters k = 0; % indicates kth model parameter for k = 1:np % initialization for j = 1:ns modelparameters(k,j) = (Pmax(k)-Pmin(k))*rand(1) + Pmin(k); % randomly distribute model parameters modelparameterchanges(k,j) = (Vmax(k)-Vmin(k))*rand(1) + Vmin(k); % randomly distribute change in model parameters end end for i = 2:iter for j = 1:ns x = modelparameters(:,j); % calculate cost function costfunction = 105*(x(2)-x(1)^2)^2 + (1-x(1))^2; if costfunction setbestcostfunction(j) = costfunction; end 4、退火算法 sjs %加载敌方 100 个目标的数据 d1=[70,40]; sj=[d1;sj;d1]; sj=sj*pi/180; %距离矩阵 d d=zeros(102); for i=1:101 for j=i+1:102 temp=cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2)); d(i,j)=6370*acos(temp); end end d=d+d'; S0=[];Sum=inf; rand('state',sum(clock)); for j=1:1000 S=[1 1+randperm(100),102]; temp=0; for i=1:101 temp=temp+d(S(i),S(i+1)); end if temp end end e=0.1^30;L=20000;at=0.999;T=1; %退火过程 for k=1:L %产生新解 c=2+floor(100*rand(1,2)); c=sort(c); c1=c(1);c2=c(2); %计算代价函数值 df=d(S0(c1-1),S0(c2))+d(S0(c1),S0(c2+1))-d(S0(c1-1),S0(c1))-d(S0(c2),S0(c2+1)); %接受准则 if df<0 S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)]; Sum=Sum+df; elseif exp(-df/T)>rand(1) S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)]; Sum=Sum+df; end T=T*at; if T end end % 输出巡航路径及路径长度 S0,Sum path=S0; xx=sj(path,1);yy=sj(path,2); plot(xx,yy,'-o')下载本文
