目录
元胞数组：相当于c中的结构体、c++中的对象
ACO（蚁群算法）
旅行商问题(TSP)优化
ACATSP.m
function [R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=ACATSP(City,Ite,Ant_num,Alpha,Beta,Rho,Q)
%% 主要符号说明
% City n个城市的坐标，n×2的矩阵
% Ite 最大迭代次数
% Ant_num 蚂蚁个数
% Alpha 表征信息素重要程度的参数
% Beta 表征启发式因子重要程度的参数
% Eta 启发因子，这里设为距离的倒数
% Tau 信息素矩阵
% Rho 信息素蒸发系数
% Q 信息素增加强度系数
% R_best 各代最佳路线
% L_best 各代最佳路线的长度
%% 第一步：变量初始化
City_num = size(City,1);%n表示问题的规模（城市个数）
Distance = zeros(City_num,City_num);%D表示完全图的赋权邻接矩阵
for i = 1:City_num
for j = 1:City_num
if i ~= j
Distance(i,j) = sqrt( (City(i,1)-City(j,1))^2 + (City(i,2)-City(j,2))^2 );% 计算任意两点间距离
else
Distance(i,j) = eps; % i=j时不计算，应该为0，但后面的启发因子要取倒数，用eps（浮点相对精度）表示
end
Distance(j,i) = Distance(i,j); % 对称矩阵
end
end
Eta = 1./Distance; % Eta为启发因子，这里设为距离的倒数
Tau = ones(City_num,City_num); % Tau为信息素矩阵
Ite_num = 1; % 迭代计数器，记录迭代次数
R_rec = zeros(Ant_num,City_num); % 存储并记录路径的生成
R_best = zeros(Ite,City_num); % 各代最佳路线
L_best = inf.*ones(Ite,1); % 各代最佳路线的长度
L_ave = zeros(Ite,1); % 各代路线的平均长度
while Ite_num<=Ite % 停止条件之一：达到最大迭代次数，停止
%% 第二步：将Ant_num只蚂蚁放到City_num个城市上
Init_ant_position = []; % 将初始状态下的蚂蚁随机分配到各个城市的临时矩阵。
for i = 1:( ceil(Ant_num/City_num) )
Init_ant_position = [Init_ant_position,randperm(City_num)];
% 每次迭代将蚂蚁随机分配到所有城市。生成一个1行多列（由于ceiling向上取整，则列数大于等于蚂蚁个数）的矩阵。
end
R_rec(:,1) = (Init_ant_position(1,1:Ant_num))';
% 矩阵转置后变成 Ant_num行-1列的一个一个初始矩阵，每一行代表一只蚂蚁，每个值代表当前蚂蚁所在城市代码。
%% 第三步：Ant_num只蚂蚁按概率函数选择下一座城市，完成各自的周游
% 这里说明下：这是个两重的大循环，外层是城市，内层是蚂蚁。
% 也就是说每次迭代每只蚂蚁向前走一步，而不是一只蚂蚁走完所有城市再开始下一只。
for j = 2:City_num % 所在城市不计算
for i = 1:Ant_num % 对每一只蚂蚁
City_visited = R_rec(i,1:(j-1)); % 记录已访问的城市，避免重复访问
City_unvisited = zeros(1,(City_num-j+1)); % 待访问的城市，初始为空
P = City_unvisited; % 待访问城市的选择概率分布，我猜这里作者弄了个简便写法，其实只是想弄一个同型矩阵。
count = 1; % 待访问城市 City_unvisited 的下标计数器
% 统计未去过的城市
for k = 1:City_num
if isempty( find( City_visited == k, 1 ) )
% 如果去过k城市，则find不为空,find(x,1)的意思是找到第一个就结束，是一个提高运算性能的写法。
% 这句话是为了避免重复去同一个城市。
City_unvisited(count) = k; % 如果if判断为真，说明第k 个城市没有去过。
count = count+1; % 下标计数器加1
end
end
% 下面计算待选城市的概率分布
for k = 1:length(City_unvisited)
P(k) = ( Tau( City_visited(end), City_unvisited(k) )^Alpha )*...
( Eta( City_visited(end), City_unvisited(k) )^Beta );
end
P=P/(sum(P));
% 按概率原则选取下一个城市
P_cum = cumsum(P); % cumsum函数是一个比较特殊的求和函数，这里是得到P 的累积概率矩阵。
Select = find(P_cum>=rand); % 若计算的概率大于原来的就选择这条路线
To_visit = City_unvisited(Select(1)); % Select(1)的意思是选中第一个累积概率大于rand随机数的城市
R_rec(i,j) = To_visit; % R_rec(i,j) 是指第i只蚂蚁，第将要去j步将要去的那个城市，循环结束后得到每只蚂蚁的路径
end
end
% 如果不是第一次循环，则将最优路径赋给路径记录矩阵的第一行
if Ite_num >= 2
R_rec(1,:) = R_best(Ite_num-1,:);
end
%% 第四步：记录本次迭代最佳路线
Len = zeros(Ant_num,1); % length 距离矩阵，初始为0。记录每只蚂蚁当前路径的总距离。
for i=1:Ant_num
R_temp = R_rec(i,:); % 取得第i 只蚂蚁的路径
% 计算第i只蚂蚁走过的总距离
for j = 1:(City_num-1)
Len(i) = Len(i) + Distance( R_temp(j),R_temp(j+1) ); % 原距离加上第j个城市到第j+1个城市的距离
end
Len(i)=Len(i)+Distance(R_temp(1),R_temp(City_num)); % 一轮下来后走过的距离
end
[L_best(Ite_num), index] = min(Len); % 最佳距离取最小
R_best(Ite_num,:) = R_rec(index(1), :);
% 此轮迭代后的最佳路线。为什么是index(1)，这是严谨写法：因为min求出后如果有多个最小值则index不唯一。
L_ave(Ite_num) = mean(Len); % 此轮迭代后的平均距离
Ite_num=Ite_num+1; % 迭代继续
%% 第五步：更新信息素
Delta_Tau = zeros(City_num, City_num); % 开始时信息素为n*n的0矩阵
for i = 1:Ant_num
for j = 1:(City_num-1)
Delta_Tau(R_rec(i,j), R_rec(i,j+1)) = Delta_Tau(R_rec(i,j), R_rec(i,j+1))+Q/Len(i);
%此次循环在路径（i，j）上的信息素增量
end
Delta_Tau(R_rec(i,City_num), R_rec(i,1)) = Delta_Tau(R_rec(i,City_num), R_rec(i,1))+Q/Len(i);
%此次循环在整个路径上的信息素增量
end
Tau = (1-Rho).*Tau + Delta_Tau; %考虑信息素挥发，更新后的信息素 Rho 信息素蒸发系数
%% 第六步：禁忌表清零
R_rec=zeros(Ant_num,City_num); %%直到最大迭代次数
end
%% 第七步：输出结果
Pos = find(L_best==min(L_best)); % 找到最佳路径（非0为真）
Shortest_Route=R_best(Pos(1), :) % 最大迭代次数后最佳路径
Shortest_Length=L_best(Pos(1) ) % 最大迭代次数后最短距离
figure % 绘制第一个子图形
DrawRoute(City,Shortest_Route) % 画路线图的子函数
figure % 绘制第二个子图形
plot(L_best)
hold on % 保持图形
plot(L_ave,'r')
title('平均距离和最短距离') % 标题DrawRoute.m
function DrawRoute(C,R)
%% 画路线图的子函数
% C Coordinate 节点坐标，由一个N×2的矩阵存储
% R Route 路线
N=length(R);
scatter(C(:,1),C(:,2));
hold on
plot([C(R(1),1),C(R(N),1)],[C(R(1),2),C(R(N),2)],'g')
for ii=2:N
plot([C(R(ii-1),1),C(R(ii),1)],[C(R(ii-1),2),C(R(ii),2)],'g')
end
title('旅行商问题优化结果 ')
hold offrun.m
City = rand(20,2).*randi([1,100],20,2);
Ite = 200;
Ant_num = 40;
Alpha = 0.7;
Beta = 0.9;
Rho = 0.3;
Q = 100;
[R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=ACATSP(City,Ite,Ant_num,Alpha,Beta,Rho,Q)citys_data.mat: 我的资源
最小二乘法
L-M法
lmm.m
function [x,val,k]=lmm(Fk,JFk,x0)
% 用L-M方法求解非线性方程组:F(x)=0
maxk=100; %最大迭代次数
rho=0.55; sigma=0.4;muk=norm(feval(Fk,x0));
k=0; epsilon=1e-6; n=length(x0);
while(k<maxk)
fk=feval(Fk,x0); %计算函数值
jfk=feval(JFk,x0); %计算Jacobi矩阵
gk=jfk'*fk;
dk=-(jfk'*jfk+muk*eye(n))\gk; %计算搜索方向
if(norm(gk)<epsilon), break; end %检验终止准则
m=0;mk=0;
while(m<20) %用Armijo准则搜索求步长
newf=0.5*norm(feval(Fk,x0+rho^m*dk))^2;
oldf=0.5*norm(feval(Fk,x0))^2;
if(newf<oldf+sigma*rho^m*gk'*dk)
mk=m; break;
end
m=m+1;
end
x0=x0+rho^mk*dk;
muk=norm(feval(Fk,x0));
k=k+1;
end
x=x0;
val=0.5*muk^2;FK.m
%% 数据拟合
%y(1)=23*x(1)+12*x(2)-2*x(3)-19;
%y(2)=-2*x(1)+32*x(2)-67*x(3)-219;
%y(3)=-18*x(1)+45*x(2)-3*x(3)-29;
%% 6
t=0.5;
% t=1;
% t=5;
y(1)=x(1)^2+x(2)^2+x(3)^2-1;
y(2)=x(1)+x(2)+x(3)-1;
y(3)=x(1)^2+x(2)^2+(x(3)-2)^2-1;
y(4)=x(1)+x(2)-x(3)+1;
y(5)=x(1)^3+3*x(2)^2+(5*x(3)-x(1)+1)^2-36*t;
y=y(:);JFK.m
function JF=JFk(x)
%JF=[1-0.7*cos(x(1)),0.2*sin(x(2));
%0.7*sin(x(1)),1+0.2*cos(x(2))];
%% 数据拟合
%JF=[23,12,-2;
% -2,32,-67;
% -18,45,-3];
%% 6
JF=[2*x(1),2*x(2),2*x(3);
1,1,1
2*x(1),2*x(2),2*(x(3)-2);
1,1,-1;
3*x(1)^2,6*x(2),10*(5*x(3)-x(1)+1)];run.m
%x0=[0,0]';
%x0=[0,0,0]';
x0=[0,0,1]';
[x,val,k]=lmm('Fk','JFk',x0)高斯-牛顿法
GaussNewton.m
function [X,Y]=GaussNewton
clear;
clc;
tic
%% 高斯牛顿算法，X0为初始点，e为迭代误差值
x0=[0;0]; %初始点
e=0.001; %阈值
syms x1 x2
F(1)=x1-0.7*sin(x1)-0.2*cos(x2); %目标函数值
F(2)=x2-0.7*cos(x1)+0.2*sin(x2);
A=jacobian(F,[x1,x2]); %求目标函数的雅可比矩阵；
H=A.'*A; %求A.'*A，这在后续会使用到
H=-inv(H); %对H求逆并取相反数
%% 进行迭代循环
for k=1:20
f=double(subs(F,[x1 x2],x0'))'; %对目标函数赋k次迭代初值
h=double(subs(H,[x1 x2],x0')); %对H赋k次迭代初值
a=double(subs(A.',[x1 x2],x0')); %对a,'赋k次迭代初值
x=x0+h*a*f; %求解k次迭代的解，即使目标函数最优的解
delet=sqrt(sum((x-x0).^2)); %误差值
if delet<e %判断误差是否符合阈值
fprintf('迭代次数为：%d\n',k)
break
else
x0=x; %不符合阈值，将K次迭代求解的值赋给K+1次迭代的初值
end
end
X=x; %返回最优解
Y=double(subs(F,[x1 x2],x'))'; %返回最优值
toc
end外罚函数+粒子群
PSO.m
function [gbest,gbestfitness,yy]=PSO(vmin,vmax,popmin,popmax,popsize,D,maxgen)
c1 = 1.49445;
c2 = 1.49445;
lu=[popmin*ones(1,D);popmax*ones(1,D)]; %%%%把边界转成向量形式
vlu=[-vmin*ones(1,D);vmax*ones(1,D)]; %%%%把边界转成向量形式
v=initpop(vlu,popsize,D); %%%%初始化速度
pop=initpop(lu,popsize,D); %%%%初始化种群
fitness=fun(pop); %%%%计算函数值
pbest=pop; %%%%初始化pbest
pbestfitness=fitness; %%%%初始化pbest函数值
%%%%进入循环
for i=1:maxgen
[gbestfitness,gbestindex]=min(fitness); %%%%%找到gbest函数值
gbest=pop(gbestindex,:); %%%%%gbest对应的个体
for j=1:popsize
v(j,:)=v(j,:)+c1*rand*(pbest(j,:)-pop(j,:))+c2*rand*(gbest-pop(j,:)); %%%%速度更新
v(j,find(v(j,:)>vmax))=vmax; %%%%边界处理，防止越界
v(j,find(v(j,:)<vmin))=vmin;
newpop(j,:)=pop(j,:)+v(j,:); %%%%%%位置更新
newpop(j,find(newpop(j,:)>popmax))=popmax; %%%%边界处理，防止越界
newpop(j,find(newpop(j,:)<popmin))=popmin;
newfitness(j)=fun(newpop(j,:)); %%%%%%计算新个体人函数值
%%更新过程
if newfitness(j)<fitness(j) %%%%更新个体%%新的函数值和原来的函数值进行对比
pop(j,:)=newpop(j,:);
fitness(j)=newfitness(j);
end
if newfitness(j)<pbestfitness(j) %%%%更新pbest%%新的函数值和pbest进行对比
pbest(j,:)=newpop(j,:);
pbestfitness(j)=newfitness(j);
end
end
yy(i)=gbestfitness; %%%%%%记录每一次循环的函数值，可以用于画图
endpenalty.m
function phvalue=penalty(x)
phvalue=0;
delta=10^15; %%%惩罚因子
%%%所有的不等式约束--》先变成大于等于0的形式
%g(1)=1-2*x(1)^2-x(2)^2;
%% 1.(3)
%g(1)=2-2*x(1)-x(2);
%g(2)=x(2)-1;
%% 1.(4)
%g(1)=x(1)+x(2)^2-1;
%g(2)=x(1)+x(2);
%% 2.(1)
%g(1)=x(1)^2-x(2);
%g(2)=x(1);
%% 2.(2)
%g(1)=-2*x(1)-x(2)+2;
%g(2)=x(2)-1;
%% 2.(3)
%g(1)=1-2*x(1)^2-x(2)^2;
%% 2.(4)
%g(1)=x;
%% 3.(1)
%g(1)=x(1);
%% 3.(2)
%g(1)=0;
%% 3.(3)
%g(1)=x(1);
%g(2)=x(2)-1;
%% 3.(4)
%g(1)=(x(1)-2)^2-x(2);
%% 4
%g(1)=-2*x(1)+x(2)+3;
%% 5
%g=[];
%% 6.(1)
g(1)=-0.25*x(1)^2-x(2)^2+1
%% 6.(2)
%g(1)=x(1);
%g(2)=x(2);
%g(3)=x(3);
%不等式约束+惩罚值
for k=1:length(g)
phvalue=phvalue+ delta*g(k)^2*getH(g(k));
end
%%所有的等式约束，如果没有就是即为空
%geq=[];
%% 1.(1)
%geq=x(1)^2+x(2)^2-1; % 把常数项移到左边，等式右边=0
%% 1.(2)
%geq=x(1)+x(2)-1;
%% 1.(3)/ 2.(1)(2)(3)(4)/3.(3)
%geq=[];
%% 3.(2)
%geq=x(1)+x(2)-1;
%% 3.(4)
%geq=2*x(1)-x(2)-1;
%% 5
%geq=x(1)+x(2)-1;
%% 6.(1)
geq= x(1)-2*x(2)+1;
%% 6.(2)
%geq(1)=8*x(1)+14*x(2)+7*x(3)-56;
%geq(2)=x(1)^2+x(2)^2+x(3)^2-25;
% 等式约束
for k=1:length(geq)
phvalue=phvalue+delta*geq(k)^2*getHeq(geq(k));
end
initpop.m
function pop=initpop(lu,N,D)
for i=1:N
for j=1:D
pop(i,j)=lu(1,D)+rand*(lu(2,D)-lu(1,D));
end
endgetH.m
function H=getH(g)
if g>=0 %%%%如果满足不等式约束条件，不进行约束，定义为0，否则进行约束定义为1
H=0;
else
H=1;
endgetHeq.m
function H=getHeq(geq)
if geq==0 %%%%如果满足等式约束条件，不进行约束，定义为0，否则进行约束定义为1
H=0;
else
H=1;
endfun.m
function y = fun(x)
%函数用于计算粒子适应度值
%x input 输入粒子
%y output 粒子适应度值
for i=1:size(x,1)
%y(i)=2*x(i,1)+3*x(i,2)+penalty(x(i,:));
%% 1.(1)
%y(i)=-x(1)-x(2)+penalty(x(i,:));
%% 1.（2）
%y(i)=x(1)^2+x(2)^2+penalty(x(i,:));
%% 1.(3)
%y(i)=x(1)^2+x(2)^2+penalty(x(i,:));
%% 1.(4)
%y(i)=-x(1)*x(2)+penalty(x(i,:));
%% 2.(1)
%y(i)=x(1)+x(2)+penalty(x(i,:));
%% 2.(2)
%y(i)=x(1)^2+x(2)+penalty(x(i,:));
%% 2.(3)
%y(i)=2*x(1)+3*x(2)+penalty(x(i,:));
%% 2.（4）
%y=(x+1).^2+penalty(x(i,:))
%% 3.(1)
%y(i)=x(1)^2+x(2)^2+penalty(x(i,:));
%% 3.(2)
%y(i)=1/2*x(1)^2+1/6*x(2)^2+penalty(x(i,:));
%% 3.(3)
%y(i)=x(1)+1/3*(x(2)+1)^2+penalty(x(i,:));
%% 3.(4)
%y(i)=(x(1)-2)^2+(x(2)-3)^2+penalty(x(i,:));
%% 4
%y(i)=x(1)*x(2)+penalty(x(i,:));
%% 5
%y(i)=x(1)^3+x(2)^3+penalty(x(i,:));
%% 6.（1）
y(i)=(x(1)-2)^2+(x(2)-1)^2+penalty(x(i,:));
%% 6.(2)
%y(i)=1000-x(1)^2-2*x(2)^2-x(3)^2-x(1)*x(2)-x(1)*x(3)+penalty(x(i,:));
endrun.m
clear
clc
maxgen=100; %%%%最大循环次数
popsize=100; %%%%种群规模
%D=1; %2.（4）
%D=3; %6.(2)
D=2; %%%问题维数，即变量个数
vmin=-0.5; %%%飞行速度上下界
vmax=0.5;
popmax=100; %%%可行域，搜索区间上下界
popmin=-100;
[gbest,gbestfitness,yy]=PSO(vmin,vmax,popmin,popmax,popsize,D,maxgen);
plot(yy)线性规划问题-单纯形法
linp.m
function [x,f,it]=linp(A,b,c) %输出x为最优解，f为最优值，it为迭代次数
b=b(:); %变为列向量: b=b';
it=0; %迭代次数
[m,n]=size(A); %%%m=3,n=2
x=zeros(1,n+length(b)); %x初始为0向量（最终存储最优值 ）
A=[A eye(length(b)) b]; %化为标准型，A b合一块
c=[c zeros(1,length(b)+1)]; %同上
while ~all(c(1:length(c)-1)>=0) %并非所有的c中前length(c)-1个元素都大于等于0时进入循环
d=find(c<0); %d(1)----第一个负数元素列坐标
e=find(A(:,d(1))>0); %e包含的d(1)列中正元素的行坐标
g=[];
for ii=1:length(e)
A(e(ii),n+length(b)+1);
A(e(ii),d(1));
g=[g A(e(ii),n+length(b)+1)/A(e(ii),d(1))];
end
h=find(g==min(g)); %选择离基变量
p=A(e(h),d(1));
for ii=1:n+length(b)+1
A(e(h),ii)=A(e(h),ii)/p; %离基变量 A(e(h),d(1)),对该行进行操作
end
j=-c(d(1))/A(e(h),d(1));
for ii=1:n+length(b)+1
c(ii)=j*A(e(h),ii)+c(ii); %%%%对c操作
end
for ii=[1:e(h)-1,e(h)+1:m]
j=-A(ii,d(1))/A(e(h),d(1));
for kk=1:n+length(b)+1
A(ii,kk)=j*A(e(h),kk)+A(ii,kk);
end
end %%%%截止，对A的操作完成
it=it+1;
end
o=[];
for ii=1:n
if all(A(:,ii)>=0)&&sum(A(:,ii))==1
o=[o ii];
end %x解的列坐标
end
for ii=1:length(o)
for kk=1:m %x解的列坐标
if A(kk,o(ii))==1
x(o(ii))=A(kk,n+length(b)+1);%对x解进行整理
end
end
end
x=x(:);
f=-c(n+length(b)+1);
end
main.m
A=[-1 1;1 2;3 1];
b=[2 10 15];
c=[-2 -3];
[x,f,it]=linp(A,b,c)BFGS+乘子法
multphr.m
function [x,mu,lambda,output]=multphr(fun,hf,gf,dfun,dhf,dgf,x0)
% 功能: 用乘子法解一般约束问题: min f(x), s.t. h(x)=0, g(x)>=0
%输入: x0是初始点, fun, dfun分别是目标函数及其梯度；
% hf, dhf分别是等式约束（向量）函数及其Jacobi矩阵的转置；
% gf, dgf分别是不等式约束（向量）函数及其Jacobi矩阵的转置；
%输出: x是近似最优点，mu, lambda分别是相应于等式约束和不
% 等式约束的乘子向量；output是结构变量，输出近似极小值f, 迭
% 代次数
maxk=500; %最大迭代次数
sigma=2.0; %罚因子
eta=2.0; theta=0.8; %PHR算法中的实参数
k=0; ink=0; %k, ink分别是外迭代和内迭代次数
epsilon=1e-5; %终止误差值
x=x0; he=feval(hf,x); gi=feval(gf,x);
n=length(x); l=length(he); m=length(gi);
%选取乘子向量的初始值
mu=0.1*ones(l,1); lambda=0.1*ones(m,1);
btak=10; btaold=10; %用来检验终止条件的两个值
while(btak>epsilon & k<maxk)
%调用BFGS算法程序求解无约束子问题
[x,ival,ik]=bfgs('mpsi','dmpsi',x0,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma);
ink=ink+ik;
he=feval(hf,x); gi=feval(gf,x);
btak=0.0;
for (i=1:l), btak=btak+he(i)^2; end
for i=1:m
temp=min(gi(i),lambda(i)/sigma);
btak=btak+temp^2;
end
btak=sqrt(btak);
if btak>epsilon
if(k>=2 & btak > theta*btaold)
sigma=eta*sigma;
end
%更新乘子向量
for (i=1:l), mu(i)=mu(i)-sigma*he(i); end
for (i=1:m)
lambda(i)=max(0.0,lambda(i)-sigma*gi(i));
end
end
k=k+1;
btaold=btak;
x0=x;
end
f=feval(fun,x);
output.fval=f;
output.iter=k;
output.inner_iter=ink;
output.bta=btak;
bfgs.m
function [x,val,k]=bfgs(fun,gfun,x0,varargin)
%功能: 用BFGS算法求解无约束问题: min f(x)
%输入: x0是初始点, fun, gfun分别是目标函数及其梯度;
% varargin是输入的可变参数变量, 简单调用bfgs时可以忽略它,
% 但若其它程序循环调用该程序时将发挥重要的作用
%输出: x, val分别是近似最优点和最优值, k是迭代次数.
maxk=500; %给出最大迭代次数
rho=0.55; sigma1=0.4; epsilon1=1e-5;
k=0; n=length(x0);
Bk=eye(n); %Bk=feval('Hess',x0);
while(k<maxk)
gk=feval(gfun,x0,varargin{:}); %计算梯度
if(norm(gk)<epsilon1), break; end %检验终止准则
dk=-Bk\gk; %解方程组, 计算搜索方向
m=0; mk=0;
while(m<20) % 用Armijo搜索求步长
newf=feval(fun,x0+rho^m*dk,varargin{:});
oldf=feval(fun,x0,varargin{:});
if(newf<oldf+sigma1*rho^m*gk'*dk)
mk=m; break;
end
m=m+1;
end
%BFGS校正
x=x0+rho^mk*dk;
sk=x-x0; yk=feval(gfun,x,varargin{:})-gk;
if(yk'*sk>0)
Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);
end
k=k+1; x0=x;
end
val=feval(fun,x0,varargin{:}); mpsi.m
%%%%%%%%%%%%%%%%%%%%%%%增广拉格朗日函数%%%%%%%%%%%%%%%%%%%
function psi=mpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma)
f=feval(fun,x); he=feval(hf,x); gi=feval(gf,x);
l=length(he); m=length(gi);
psi=f; s1=0.0;
for(i=1:l)
psi=psi-he(i)*mu(i);
s1=s1+he(i)^2;
end
psi=psi+0.5*sigma*s1;
s2=0.0;
for(i=1:m)
s3=max(0.0, lambda(i) - sigma*gi(i));
s2=s2+s3^2-lambda(i)^2;
end
psi=psi+s2/(2.0*sigma);
dmpsi.m
%%%%%%%%%%%%%%%%%%%%%%%增广拉格朗日函数的梯度%%%%%%%%%%%%%%%%%%%
function dpsi=dmpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma)
dpsi=feval(dfun,x);
he=feval(hf,x); gi=feval(gf,x);
dhe=feval(dhf,x); dgi=feval(dgf,x);
l=length(he); m=length(gi);
for(i=1:l)
dpsi=dpsi+(sigma*he(i)-mu(i))*dhe(:,i);
end
for(i=1:m)
dpsi=dpsi+(sigma*gi(i)-lambda(i))*dgi(:,i);
end
f1.m
function f=f1(x)
%% 例题和习题9. 6(1)
%f=(x(1)-2.0)^2+(x(2)-1.0)^2;
%% 9.6(2)
f=1000-x(1)^2-2*x(2)^2-x(3)^2-x(1)*x(2)-x(1)*x(3);df1.m
function g=df1(x)
%% 例题和习题9. 6(1)
%g=[2.0*(x(1)-2.0), 2.0*(x(2)-1.0)]';
%% 9.6(2)
g=[-2*x(1)-x(2)-x(3),-4*x(2)-x(1),-2*x(3)-x(1)]';g1.m
%不等式约束 xxx>=0
function gi=g1(x)
%% 例题 和 9.6（1）
%gi=-0.25*x(1)^2-x(2)^2+1;
%% 习题9.6(2)
gi(1)=x(1);
gi(2)=x(2);
gi(3)=x(3);dg1.m
function dgi=dg1(x)
%% 例题 和 9.6（1）
%dgi=[-0.5*x(1),-2.0*x(2)]';
%% 9.6(2)
dgi=[1,0,0;
0,1,0;
0,0,1]'h1.m
% 等式约束 xxx=0
function he=h1(x)
%% 例题 和 9.6（1）
%he=x(1)-2.0*x(2)+1.0;
%% 习题9.6(2)
he(1)=8*x(1)+14*x(2)+7*x(3)-56;
he(2)=x(1)^2+x(2)^2+x(3)^2-25;dh1.m
function dhe=dh1(x)
%% 例题 和 9.6（1）
%dhe=[1.0,-2.0]';
%% 9.6（2）
dhe=[8.0,14.0,7.0;
2*x(1),2*x(2),2*x(3);
0,0,0]'run.m
%% 例题
%x0=[3,3]';
%% 9.6（1）
%x0=[2,2]';
%% 9.6(2)
x0=[2,2,2]';
[x,mu,lambda,output]=multphr('f1','h1','g1','df1','dh1','dg1',x0)二次规划
拉格朗日方法
qlag.m
function [x,lam,fval]=qlag(H,A,b,c)
% 功能：用拉格朗日方法求解等式约束二次规划
% min f(x)=0.5*x'Hx+c'x, s.t.Ax=b
% 输入：H，c分别是目标函数的矩阵和向量，A,b分别是
% 约束条件中的矩阵和向量
% 输出：(x,lam)是 KT 点，fval 是最优值
IH=inv(H);
AHA=A*IH*A';
IAHA=inv(AHA);
AIH=A*IH;
G=IH-AIH'*IAHA*AIH;
B=IAHA*AIH;
C=-IAHA;
x=B'*b-G*c;
lam=B*c-C*b;
fval=0.5*x'*H*x+c'*x
run.m
H=[3 -2 0;-2 2 -2;0 -2 1];
c=[1 1 1]';
A=[1 2 1];
b=[4]';
[x,lam]=qlag(H,A,b,c)quadprog命令
[x,fval,exitflag,output,lambda]=quadprog(c,A,b,Aeq,beq,lb,ub,x0,options)
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
