目录

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_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 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;
        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('迭代次数为：%dn',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;                             %%%%%%记录每一次循环的函数值，可以用于画图
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:N
    for j=1:D
        pop(i,j)=lu(1,D)+rand*(lu(2,D)-lu(1,D));
    end
end

getH.m

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

getHeq.m

function H=getHeq(geq)
if geq==0       %%%%如果满足等式约束条件，不进行约束，定义为0，否则进行约束定义为1
   H=0;
else
   H=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
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=-Bkgk;  %解方程组, 计算搜索方向
    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)

最后

以上就是端庄万宝路为你收集整理的约束问题的最优化ACO（蚁群算法）最小二乘法外罚函数+粒子群线性规划问题-单纯形法BFGS+乘子法二次规划的全部内容，希望文章能够帮你解决约束问题的最优化ACO（蚁群算法）最小二乘法外罚函数+粒子群线性规划问题-单纯形法BFGS+乘子法二次规划所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。