部分MATLAB函数实现

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 genFitnV=ranking(-ObjV); %分配适应度函数值

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 ccpoint=ceil(rand*(2*L-1)); %取得一个1到2L-1的整数/ceil表示+1法取整，找到交叉的位置

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)mm(N,:)=zeros(1,2*L); %最后一行不变异，强制赋0

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))flag=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 　　bestmodelparameters(:,j) = modelparameters(:,j);

　　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 tempS0=S;Sum=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 Tbreak;

end

end

% 输出巡航路径及路径长度

S0,Sum

path=S0;

xx=sj(path,1);yy=sj(path,2);

plot(xx,yy,'-o')