阳山网站建设,《网站建设与管理》方案,网站集约化建设困难,wordpress nextgen gallery1. 理论部分
2. 序列二次规划算法代码及解析
3.完整代码
1.理论部分
a.约束优化问题的极值条件
库恩塔克条件(Kuhn-Tucker conditions#xff0c;KT条件)是确定某点为极值点的必要条件。如果所讨论的规划是凸规划#xff0c;那么库恩-塔克条件也是充分条件。
#xff…1. 理论部分
2. 序列二次规划算法代码及解析
3.完整代码
1.理论部分
a.约束优化问题的极值条件
库恩塔克条件(Kuhn-Tucker conditionsKT条件)是确定某点为极值点的必要条件。如果所讨论的规划是凸规划那么库恩-塔克条件也是充分条件。
1判断数值x处起作用的约束计算x处各约束的取值是否为0找出起作用的约束条件。并非所有约束条件都起作用即并非落在所有约束条件的边界交集处。对于起作用的边界条件满足边界函数g(x)或h(x)为0
2数值x处目标函数导数f(x)和起作用的约束导数g(x)的数值
3代入KT条件ex. 求解拉格朗日乘子λ的大小判断λ是否大于0若大于零则为极值点。
4若不满足条件继续迭代。
具体做法如下 其中λ为拉格朗日乘子当某一点x使λ 0存在时称满足KT条件则该点x是极值点。 约束优化问题的应用示例如下图 b.约束优化方法及迭代
此处介绍和使用外惩罚函数法/外点法
外点法是一种从随机一点一般在可行域外出发逐渐移动到可行区域的方法。 对于约束优化问题 构造外惩罚函数 即当g(x)不满足条件时该项有大于0的值h(x)不满足条件时该项也有值加和的越大代表该数值点距离距离约束的可行域越远惩罚项的值越大惩罚作用越明显。。其中Mk为外罚因子M 0 M 1 M 2 . . . M k 为递增的正数序列一般M01递增系数c5~8。
迭代终止满足两个判断条件中一个即可
1. Q为当前数值点多个约束函数中最大的值即满足了所有的约束函数条件。
2. 且 R为外罚因子控制量M超出外罚因子的极限R
外点法的迭代图如下 外点法的示例如下图 a. 定义变量
vars 2; % 自变量数量
cons 1; % 约束数量
maxIter100; % 最大迭代次数
x [-1;4]; % 初始猜测值
l 0; % 拉格朗日乘子
H eye(vars,vars); % Hessian 矩阵假设为I矩阵
目标函数为 及目标函数的导数为 function y f(x)y x(1)^4 - 2*x(2)*x(1)^2 x(2)^2 x(1)^2 - 2*x(1)5;
endfunction y gradf(x)y(1) 2*x(1)-4*x(1)*x(2)4*x(1)^3-2;y(2) -2*x(1)^2 2*x(2);
end
约束项不等式为 以及其导数为 function y g(x)y(1) -(x(1)0.25)^20.75*x(2);
endfunction y gradg(x)y(1,1) -2*x(1)-1/2;y(1,2) 3/4;
end
求解KKT条件
function y SolveKKT(gradfEval,gradgEval,gEval,Hessian)A [Hessian -gradgEval;gradgEval 0]; %对称矩阵约束函数的梯度矩阵b [-gradfEval -gEval];y A\b; %方程 Ay b 解出y值分别对应最优的x值的梯度和最优的拉格朗日乘子的值%即此时Ay-b 0kkt条件成立
end
判断数值点x是否满足约束条件
function [gViol,lViol] Viols(gEval,l)%约束条件函数值拉格朗日乘子lgViol[];lViol[];for i 1:length(gEval)if gEval(i)0 %如果约束条件函数值小于零即不满足约束条件gViol(i)gEval(i); %将约束条件函数值和拉格朗日乘子l放入列表lViol(i)l(i);end end
end
处理不满足约束的情况构造惩罚函数lViol对应Mk累加得到惩罚函数
function y Penalty(fEval,gViol,lViol)
%代入该点处的函数值不满足约束的约束函数值以及拉格朗日乘子lsum 0;y fEval;for i 1:length(gViol) %找出不满足约束条件的约束函数值sum sum lViol(i)*abs(gViol(i)); %约束函数值g(x)求反再乘一个系数即||l_1*g(x1)l2*g(x2) ||endy y sum;
end
在初始点首先进行初始计算找出不满足约束的约束计算其损失函数。
% EVALUATE AT STARTING POINTfEval f(x);
gEval g(x);
[gViol,lViol] Viols(gEval,l);
gradfEval gradf(x);
gradgEval gradg(x);
P Penalty(fEval,gViol,lViol);
开始进入SQP算法主体,求解kkt条件得到xSol为函数x的解lSol为拉格朗日乘子的解。
for i1:maxIter %开始循环迭代% SOLVE KKT CONDITIONS FOR THE OPTIMUM OF THE QUADRATIC APPROXIMATION%求解二次优化的KKT条件sol SolveKKT(gradfEval,gradgEval,gEval,H);%代入导数f(x)g(x),约束函数值g(x)Hessian矩阵xSol sol(1:vars);lSol sol(vars1:varscons);
如果拉格朗日乘子有负的即不满足约束条件则将其设为0并重新求解方程 for j 1:length(lSol)if lSol(j)0 %如果拉格朗日乘子有一个是负的sol H\(-gradfEval); %将结果重新求解H*x -gradfEvalxSol sol(1:vars);lSol(j)0; %该约束的拉格朗日乘子置为0endend
计算新的数值点情况 xNew x xSol; %更新得到新的数值点lNew lSol; % 得到新的拉格朗日乘子fEvalNew f(xNew);%计算新点目标函数的数值gEvalNew g(xNew);%计算新点的约束函数数值gradfEvalNew gradf(xNew);%计算该点目标函数梯度gradgEvalNew gradg(xNew);%计算该点约束函数梯度[gViolNew,lViolNew] Viols(gEvalNew,lNew);%找出不满足的约束条件PNew Penalty(fEvalNew,gViolNew,lViolNew);%求解惩罚函数
如果惩罚函数增加减小步长到原来的一半 % IF PENALTY FUNCTION INCREASED, LOWER THE STEP BY 0.5while PNew-P1e-4xSol 0.5*xSol;xNew x xSol;fEvalNew f(xNew);gEvalNew g(xNew);gradfEvalNew gradf(xNew);gradgEvalNew gradg(xNew);[gViolNew,lViolNew] Viols(gEvalNew,lNew);PNew Penalty(fEvalNew,gViolNew,lViolNew);end
如果x数值变化很小停止迭代找到最优结束迭代 % STOPPING CRITERIONif norm(xNew(1:vars)-x(1:vars))1e-2breakend
更新Hessian矩阵Q为Mk1- Mk的值 % UPDATE THE HESSIANgradLEval gradLagr(gradfEval,gradgEval,lNew,vars); %拉格朗日函数的梯度 lnew not l!!!gradLEvalNew gradLagr(gradfEvalNew,gradgEvalNew,lNew,vars);%新的数值点处拉格朗日函数的梯度Q gradLEvalNew-gradLEval;%拉格朗日函数数值的差分L(x_k1) - L(x_k)dx xNew-x;%自变量x的差分HNew UpdateH(H,dx,Q);%求解H矩阵
计算拉格朗日乘子梯度的函数
function y gradLagr(gradfEval,gradgEval,l,n)y gradfEval;sum zeros(n,1);for i 1:length(l)sum sum -l(i)*gradgEval(i:n);endy y sum; %x处目标函数的梯度减去约束函数的梯度。即拉格朗日函数的梯度
end
更新Hessian矩阵gamma为Mk1- Mk的值
function y UpdateH(H,dx,gamma)%gamma是拉格朗日函数的两次差分term1(gamma*gamma) / (gamma*dx);term2 ((H*dx)*(dx*H)) / (dx*(H*dx));y H term1-term2; %更新得到新的H矩阵
end
更新所有的变量用于下次的迭代包括Hessian矩阵H数值点x惩罚函数值P目标函数值fEval目标函数值梯度gradEval约束函数值gEval约束函数值梯度gradgEval % UPDATE ALL VALUES FOR NEXT ITERATIONH HNew;fEval fEvalNew;gEval gEvalNew;gradfEval gradfEvalNew;gradgEval gradgEvalNew;P PNew;x xNew; 3.可以运行的完整代码如下
% EXAMPLE OF SQP ALGORITHMclear variables; close all; clc;
fprintf(---------------------------------------------------------------\n)
fprintf(An implementation of Sequential Quadratic Programming method\nin a nonlinear constrained optimization problem\n)
fprintf(---------------------------------------------------------------\n)% INITIAL VALUES - INPUTvars 2; % number of variables
cons 1; % number of constraints
maxIter100; % max iterations
x [-1;4]; % initial guess point
l 0; % LagrangeMultipliers vector
H eye(vars,vars); % Hessian matrix assumed to be identity% EVALUATE AT STARTING POINTfEval f(x);
gEval g(x);
[gViol,lViol] Viols(gEval,l);
gradfEval gradf(x);
gradgEval gradg(x);
P Penalty(fEval,gViol,lViol);% SQP ALGORITHMfor i1:maxIter% SOLVE KKT CONDITIONS FOR THE OPTIMUM OF THE QUADRATIC APPROXIMATIONsol SolveKKT(gradfEval,gradgEval,gEval,H);xSol sol(1:vars);lSol sol(vars1:varscons);% IF THE LAGRANGE MULTIPLIER IS NEGATIVE SET IT TO ZEROfor j 1:length(lSol)if lSol(j)0 sol H\(-gradfEval);xSol sol(1:vars);lSol(j)0;endend% EVALUATE AT NEW CANDIDATE POINTxNew x xSol; lNew lSol;fEvalNew f(xNew);gEvalNew g(xNew);gradfEvalNew gradf(xNew);gradgEvalNew gradg(xNew);[gViolNew,lViolNew] Viols(gEvalNew,lNew);PNew Penalty(fEvalNew,gViolNew,lViolNew);% IF PENALTY FUNCTION INCREASED, LOWER THE STEP BY 0.5while PNew-P1e-4xSol 0.5*xSol;xNew x xSol;fEvalNew f(xNew);gEvalNew g(xNew);gradfEvalNew gradf(xNew);gradgEvalNew gradg(xNew);[gViolNew,lViolNew] Viols(gEvalNew,lNew);PNew Penalty(fEvalNew,gViolNew,lViolNew);end% STOPPING CRITERIONif norm(xNew(1:vars)-x(1:vars))1e-2breakend% UPDATE THE HESSIANgradLEval gradLagr(gradfEval,gradgEval,lNew,vars); % lnew not l!!!gradLEvalNew gradLagr(gradfEvalNew,gradgEvalNew,lNew,vars);Q gradLEvalNew-gradLEval;dx xNew-x;HNew UpdateH(H,dx,Q);% UPDATE ALL VALUES FOR NEXT ITERATIONH HNew;fEval fEvalNew;gEval gEvalNew;gradfEval gradfEvalNew;gradgEval gradgEvalNew;P PNew;x xNew;
endfprintf(SQP: Optimum point:\n x1%10.4f\n x2%10.4f\n iterations %10.0f \n, x(1), x(2), i)% FUNCTIONS NEEDEDfunction y SolveKKT(gradfEval,gradgEval,gEval,Hessian)A [Hessian -gradgEval;gradgEval 0];b [-gradfEval -gEval];y A\b;
endfunction y f(x)y x(1)^4 - 2*x(2)*x(1)^2 x(2)^2 x(1)^2 - 2*x(1)5;
endfunction y gradf(x)y(1) 2*x(1)-4*x(1)*x(2)4*x(1)^3-2;y(2) -2*x(1)^2 2*x(2);
endfunction y gradLagr(gradfEval,gradgEval,l,n)y gradfEval;sum zeros(n,1);for i 1:length(l)sum sum -l(i)*gradgEval(i:n);endy y sum;
endfunction y gradg(x)y(1,1) -2*x(1)-1/2;y(1,2) 3/4;
endfunction y g(x)y(1) -(x(1)0.25)^20.75*x(2);
endfunction [gViol,lViol] Viols(gEval,l)gViol[];lViol[];for i 1:length(gEval)if gEval(i)0gViol(i)gEval(i);lViol(i)l(i);end end
endfunction y Penalty(fEval,gViol,lViol)sum 0;y fEval;for i 1:length(gViol)sum sum lViol(i)*abs(gViol(i));endy y sum;
endfunction y UpdateH(H,dx,gamma)term1(gamma*gamma) / (gamma*dx);term2 ((H*dx)*(dx*H)) / (dx*(H*dx));y H term1-term2;
end最后输出结果为
---------------------------------------------------------------
An implementation of Sequential Quadratic Programming method
in a nonlinear constrained optimization problem
---------------------------------------------------------------
SQP: Optimum point:x1 0.4999x2 0.7498iterations 9
