约束问题的最优化

最优化"/>

约束问题的最优化

目录

ACO（蚁群算法）

旅行商问题(TSP)优化

ACATSP.m

 DrawRoute.m

run.m 

citys_data.mat: 我的资源

最小二乘法

L-M法

lmm.m

FK.m

JFK.m

run.m

高斯-牛顿法

GaussNewton.m

外罚函数+粒子群

PSO.m

penalty.m

initpop.m

getH.m

getHeq.m

fun.m

run.m

线性规划问题-单纯形法

linp.m

main.m

BFGS+乘子法

multphr.m

bfgs.m

mpsi.m

dmpsi.m

f1.m

df1.m

g1.m

dg1.m

h1.m

dh1.m

run.m

二次规划

拉格朗日方法

qlag.m

run.m

quadprog命令

---

元胞数组：相当于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_numfor j = 1:City_numif i ~= jDistance(i,j) = sqrt( (City(i,1)-City(j,1))^2 + (City(i,2)-City(j,2))^2 );% 计算任意两点间距离elseDistance(i,j) = eps;      % i=j时不计算，应该为0，但后面的启发因子要取倒数，用eps（浮点相对精度）表示endDistance(j,i) = Distance(i,j);   % 对称矩阵end
endEta = 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_numif isempty( find( City_visited == k, 1 ) ) % 如果去过k城市，则find不为空,find(x,1)的意思是找到第一个就结束，是一个提高运算性能的写法。% 这句话是为了避免重复去同一个城市。City_unvisited(count) = k; % 如果if判断为真，说明第k 个城市没有去过。count = count+1;      % 下标计数器加1endend% 下面计算待选城市的概率分布for k = 1:length(City_unvisited)P(k) = ( Tau( City_visited(end), City_unvisited(k) )^Alpha )*...( Eta( City_visited(end), City_unvisited(k) )^Beta );endP=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步将要去的那个城市，循环结束后得到每只蚂蚁的路径endend% 如果不是第一次循环，则将最优路径赋给路径记录矩阵的第一行if Ite_num >= 2R_rec(1,:) = R_best(Ite_num-1,:);end%% 第四步：记录本次迭代最佳路线Len = zeros(Ant_num,1);     % length 距离矩阵，初始为0。记录每只蚂蚁当前路径的总距离。for i=1:Ant_numR_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个城市的距离endLen(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_numfor 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）上的信息素增量endDelta_Tau(R_rec(i,City_num), R_rec(i,1)) = Delta_Tau(R_rec(i,City_num), R_rec(i,1))+Q/Len(i);%此次循环在整个路径上的信息素增量endTau = (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:Nplot([C(R(ii-1),1),C(R(ii),1)],[C(R(ii-1),2),C(R(ii),2)],'g')
end
title('旅行商问题优化结果 ')
hold off

run.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;endm=m+1;endx0=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;%% 6t=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,12*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:20f=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)breakelsex0=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:popsizev(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);endif newfitness(j)<pbestfitness(j)               %%%%更新pbest%%新的函数值和pbest进行对比pbest(j,:)=newpop(j,:);pbestfitness(j)=newfitness(j);endendyy(i)=gbestfitness;                             %%%%%%记录每一次循环的函数值，可以用于画图
end

penalty.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:Nfor j=1:Dpop(i,j)=lu(1,D)+rand*(lu(2,D)-lu(1,D));end
end

getH.m

function H=getH(g)
if g>=0   %%%%如果满足不等式约束条件，不进行约束，定义为0，否则进行约束定义为1H=0;
elseH=1;
end

getHeq.m

function H=getHeq(geq)
if geq==0       %%%%如果满足等式约束条件，不进行约束，定义为0，否则进行约束定义为1H=0;
elseH=1;
end

fun.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,:));end

run.m

clear
clcmaxgen=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))];endh=find(g==min(g)); %选择离基变量p=A(e(h),d(1));for ii=1:n+length(b)+1A(e(h),ii)=A(e(h),ii)/p; %离基变量 A(e(h),d(1)),对该行进行操作endj=-c(d(1))/A(e(h),d(1));for ii=1:n+length(b)+1c(ii)=j*A(e(h),ii)+c(ii); %%%%对c操作endfor 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)+1A(ii,kk)=j*A(e(h),kk)+A(ii,kk);endend %%%%截止，对A的操作完成it=it+1;
end
o=[];
for ii=1:nif all(A(:,ii)>=0)&&sum(A(:,ii))==1o=[o ii];end %x解的列坐标
end
for ii=1:length(o)for kk=1:m %x解的列坐标if A(kk,o(ii))==1x(o(ii))=A(kk,n+length(b)+1);%对x解进行整理endend
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;   endfor i=1:mtemp=min(gi(i),lambda(i)/sigma);btak=btak+temp^2;endbtak=sqrt(btak);   if btak>epsilonif(k>=2 & btak > theta*btaold)sigma=eta*sigma;end%更新乘子向量for (i=1:l),  mu(i)=mu(i)-sigma*he(i);  endfor (i=1:m)lambda(i)=max(0.0,lambda(i)-sigma*gi(i));endendk=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;endm=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);endk=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)