顺序二次规划 Matlab 实现
Sequential Quadratic Programming Matlab Implementation
我的 MATLAB 代码有问题，该代码使用 SQP 算法[=21] 求解 非线性 二次问题=] （顺序二次规划） 但在 "QP-SUB PROBLEM" 部分我已经分析制定了 "num2str"error 出现，老实说，我不知道如何解决这个问题，还必须告诉你这个方法使用 KT 条件 以获得更好的解决方案.
在代码的每一节中，我都写了一条注释以便更好地理解，并且可以在下面的代码中找到带有约束的功能：
% Maximize f(x1,x2) = x1^4 -2x1^2x2 +x1^2 +x1x2^2 -2x1 +4
%
% h1(x1,x2) = x1^2 + x2^2 -2 = 0
% g1(x1,x2) = 0.25x1^2 +0.75x2^2 -1 <=0
% 0 <= x1 <= 4; 0 <= x2 <= 4
%
%--------------------------------------------------------
% The KT conditions for the QP subproblem
% are applied analytically
% There are two cases for a single inequality constraint
% Case (a) : beta = 0 g < 0
% Case (b) : beta ~= 0, g = 0
% The best solution is used
% --------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% management functions
clear % clear all variable/information in the workspace - use CAUTION
clear global % again use caution - clears global information
clc % position the cursor at the top of the screen
close % closes the figure window
format compact % avoid skipping a line when writing to the command window
warning off %#ok<WNOFF> % don't report any warnings like divide by zero etc.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% DATA ---------------------
%%% starting design
xb(1) = 3; xb(2) = 2;
it = 10; % number of iterations
%%% plot range for delx1 : -3 , +3
dx1L = -3; dx1U = +3;
dx2L = -3; dx2U = +3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Define functions
syms x1 x2 f g h
syms gradf1 gradf2 gradh1 gradh2 gradg1 gradg2
f = x1^4 - 2*x1*x1*x2 + x1*x1 + x1*x2*x2 - 2*x1 + 4;
h = x1*x1 + x2*x2 - 2;
g = 0.25*x1*x1 +0.75*x2*x2 -1;
%%% the gradient functions
gradf1 = diff(f,x1);
gradf2 = diff(f,x2);
% the hessian
hess = [diff(gradf1,x1), diff(gradf1,x2); diff(gradf2,x1), diff(gradf2,x2)];
% gradient of the constraints
gradh1 = diff(h,x1);
gradh2 = diff(h,x2);
gradg1 = diff(g,x1);
gradg2 = diff(g,x2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% graphical/symbolic solution for SLP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('***********************')
fprintf('\nSQP - Example 7.1')
fprintf('\n*********************\n')
for i = 1:it
%figure;
f1 = subs(f,{x1,x2},{xb(1),xb(2)});
g1 = subs(g,{x1,x2},{xb(1),xb(2)});
h1 = subs(h,{x1,x2},{xb(1),xb(2)});
%%% Print Information
fprintf('\n******************************')
fprintf('\nIteration : '),disp(i)
fprintf('******************************\n')
fprintf('Linearized about [x1, x2] : '),disp([xb(1) xb(2)])
fprintf('Objective function value f(x1,x2) : '),disp(f1);
fprintf('Equality constraint value value h(x1,x2) : '),disp(h1);
fprintf('Inequality constraint value value g(x1,x2) : '),disp(g1);
%fprintf('\nsolution for [delx1 delx2] : '),disp(sol')
% hold on
% calculate the value of the gradients
% f1 = subs(f,{x1,x2},{xb(1),xb(2)});
% g1 = subs(g,{x1,x2},{xb(1),xb(2)});
% h1 = subs(h,{x1,x2},{xb(1),xb(2)});
fprintf('\n-----------------------')
fprintf('\nQP - SUB PROBLEM')
fprintf('\n---------------------\n')
gf1 = double(subs(gradf1,{x1,x2},{xb(1),xb(2)}));
gf2 = double(subs(gradf2,{x1,x2},{xb(1),xb(2)}));
hess1 = double(subs(hess,{x1,x2},{xb(1),xb(2)}));
gh1 = double(subs(gradh1,{x1,x2},{xb(1),xb(2)}));
gh2 = double(subs(gradh2,{x1,x2},{xb(1),xb(2)}));
gg1 = double(subs(gradg1,{x1,x2},{xb(1),xb(2)}));
gg2 = double(subs(gradg2,{x1,x2},{xb(1),xb(2)}));
% the QP subproblem
syms dx1 dx2 % change in design
fquad = f1 + [gf1 gf2]*[dx1; dx2] + 0.5*[dx1 dx2]*hess1*[dx1 ;dx2];
hlin = h1 + [gh1 gh2]*[dx1; dx2];
glin = g1 + [gg1 gg2]*[dx1; dx2];
Fquadstr = strcat(num2str(f1),' + ',num2str(gf1), ...
'*','dx1',' + ',num2str(gf2),' * ','dx2', ...
' + 0.5*',num2str(hess1(1,1)),' * dx1^2', ...
' +',num2str(hess1(1,2)),' * dx1*dx2', ...
' + 0.5*',num2str(hess1(2,2)),' * dx2^2');
hlinstr = strcat(num2str(h1),' + ',num2str(gh1), ...
'*','dx1',' + ',num2str(gh2),' * ','dx2');
glinstr = strcat(num2str(g1),' + ',num2str(gg1), ...
'*','dx1',' + ',num2str(gg2),' * ','dx2');
fprintf('Quadratic Objective function f(x1,x2): \n'),disp(Fquadstr);
fprintf('Linearized equality h(x1,x2): '),disp(hlinstr);
fprintf('Linearized inequality g(x1,x2): '),disp(glinstr);
fprintf('\n')
% define Lagrangian for the QP problem
syms lamda beta
F = fquad + lamda*hlin + beta*glin;
fprintf('Case a: beta = 0\n');
Fnobeta = fquad + lamda*hlin;
%%% initialize best solution
dx1best = 0;
dx2best = 0;
Fbbest = 0;
%%%%%%%%%%%%%%%%%%%%%%%
%%% solve case (a) %%%
%%%%%%%%%%%%%%%%%%%%%%%
xcasea = solve(diff(Fnobeta,dx1),diff(Fnobeta,dx2),hlin);
sola = [double(xcasea.dx1) double(xcasea.dx2) double(xcasea.lamda)];
dx1a = double(xcasea.dx1);
dx2a = double(xcasea.dx2);
lamdaa = double(xcasea.lamda);
hlina = double(subs(hlin,{dx1,dx2},{dx1a,dx2a}));
glina = double(subs(glin,{dx1,dx2},{dx1a,dx2a}));
Fa = double(subs(Fnobeta,{dx1,dx2,lamda},{dx1a,dx2a,lamdaa}));
%%% results for case (a)
x1a = dx1a + xb(1);
x2a = dx2a + xb(2);
fv = double(subs(f,{x1,x2},{x1a,x2a}));
hv = double(subs(h,{x1,x2},{x1a,x2a}));
gv = double(subs(g,{x1,x2},{x1a,x2a}));
fprintf('Change in design vector: '),disp([dx1a dx2a]);
fprintf('The linearized quality constraint: '),disp(hlina);
fprintf('The linearized inequality constraint: '),disp(glina);
fprintf('New design vector: '),disp([x1a x2a]);
fprintf('The objective function: '),disp(fv);
fprintf('The equality constraint: '),disp(hv);
fprintf('The inequality constraint: '),disp(gv);
if (glina <= 0)
xb(1) = xb(1) + dx1a;
xb(2) = xb(2) + dx2a;
fbest = Fa;
dx1best = dx1a;
dx2best = dx2a;
end
%%%%%%%%%%%%%%%%%%%%%%%
%%% solve case (b) %%%
%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n Case b: g = 0\n');
xcaseb = solve(diff(F,dx1),diff(F,dx2),hlin,glin);
solb = [double(xcaseb.dx1) double(xcaseb.dx2) double(xcaseb.lamda) double(xcaseb.beta)];
dx1b = double(xcaseb.dx1);
dx2b = double(xcaseb.dx2);
betab = double(xcaseb.beta);
lamdab = double(xcaseb.lamda);
hlinb = double(subs(hlin,{dx1,dx2},{dx1b,dx2b}));
glinb = double(subs(glin,{dx1,dx2},{dx1b,dx2b}));
Fb = double(subs(F,{dx1,dx2,lamda,beta},{dx1b,dx2b,lamdab,betab}));
x1b = dx1b + xb(1);
x2b = dx2b + xb(2);
fv = double(subs(f,{x1,x2},{x1b,x2b}));
hv = double(subs(h,{x1,x2},{x1b,x2b}));
gv = double(subs(g,{x1,x2},{x1b,x2b}));
fprintf('Change in design vector: '),disp([dx1b dx2b]);
fprintf('The linearized quality constraint: '),disp(hlinb);
fprintf('The linearized inequality constraint: '),disp(glinb);
fprintf('New design vector: '),disp([x1b x2b]);
fprintf('The objective function: '),disp(fv);
fprintf('The equality constraint: '),disp(hv);
fprintf('The inequality constraint: '),disp(gv);
fprintf('Multiplier beta: '),disp(betab);
fprintf('Multiplier lamda: '),disp(lamdab);
if (betab > 0) & (Fb <= fbest)
xb(1) = x1b;
xb(2) = x2b;
dx1best = dx1b;
dx2best = dx2b;
end
%%% stopping criteria
if ([dx1best dx2best]*[dx1best dx2best]') <= 1.0e-08
fprintf('\n&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
fprintf('\nStopped: Design Not Changing')
fprintf('\n&&&&&&&&&&&&&&&&&&&&&&&&&&&&&\n\n')
break;
elseif i == it
fprintf('\n&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
fprintf('\nStpped: Number of iterations at maximum')
fprintf('\n&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&\n\n')
break;
end
end
f1, g1, h1仍然是syms变量类型。- 在应用前使用函数
double() 将它们更改为 numeric 类型
函数 num2str()
- 您已经对变量应用了
double()
gf1, gf2, gh1, gh2,gg1, gg2 and hess1
就在上面，不用碰它们
这是您应该替换的部分
Fquadstr = strcat(num2str(f1),' + ',num2str(gf1), ...
'*','dx1',' + ',num2str(gf2),' * ','dx2', ...
' + 0.5*',num2str(hess1(1,1)),' * dx1^2', ...
' +',num2str(hess1(1,2)),' * dx1*dx2', ...
' + 0.5*',num2str(hess1(2,2)),' * dx2^2');
hlinstr = strcat(num2str(h1),' + ',num2str(gh1), ...
'*','dx1',' + ',num2str(gh2),' * ','dx2');
glinstr = strcat(num2str(g1),' + ',num2str(gg1), ...
'*','dx1',' + ',num2str(gg2),' * ','dx2');
据此
% apply double() to f1
Fquadstr = strcat(num2str(double(f1)),' + ',num2str(gf1), ...
'*','dx1',' + ',num2str(gf2),' * ','dx2', ...
' + 0.5*',num2str(hess1(1,1)),' * dx1^2', ...
' +',num2str(hess1(1,2)),' * dx1*dx2', ...
' + 0.5*',num2str(hess1(2,2)),' * dx2^2');
% apply double() to h1
hlinstr = strcat(num2str(double(h1)),' + ',num2str(gh1), ...
'*','dx1',' + ',num2str(gh2),' * ','dx2');
% apply double() to g1
glinstr = strcat(num2str(double(g1)),' + ',num2str(gg1), ...
'*','dx1',' + ',num2str(gg2),' * ','dx2');
我的 MATLAB 代码有问题，该代码使用 SQP 算法[=21] 求解 非线性 二次问题=] （顺序二次规划） 但在 "QP-SUB PROBLEM" 部分我已经分析制定了 "num2str"error 出现，老实说，我不知道如何解决这个问题，还必须告诉你这个方法使用 KT 条件 以获得更好的解决方案. 在代码的每一节中，我都写了一条注释以便更好地理解，并且可以在下面的代码中找到带有约束的功能：
% Maximize f(x1,x2) = x1^4 -2x1^2x2 +x1^2 +x1x2^2 -2x1 +4
%
% h1(x1,x2) = x1^2 + x2^2 -2 = 0
% g1(x1,x2) = 0.25x1^2 +0.75x2^2 -1 <=0
% 0 <= x1 <= 4; 0 <= x2 <= 4
%
%--------------------------------------------------------
% The KT conditions for the QP subproblem
% are applied analytically
% There are two cases for a single inequality constraint
% Case (a) : beta = 0 g < 0
% Case (b) : beta ~= 0, g = 0
% The best solution is used
% --------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% management functions
clear % clear all variable/information in the workspace - use CAUTION
clear global % again use caution - clears global information
clc % position the cursor at the top of the screen
close % closes the figure window
format compact % avoid skipping a line when writing to the command window
warning off %#ok<WNOFF> % don't report any warnings like divide by zero etc.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% DATA ---------------------
%%% starting design
xb(1) = 3; xb(2) = 2;
it = 10; % number of iterations
%%% plot range for delx1 : -3 , +3
dx1L = -3; dx1U = +3;
dx2L = -3; dx2U = +3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Define functions
syms x1 x2 f g h
syms gradf1 gradf2 gradh1 gradh2 gradg1 gradg2
f = x1^4 - 2*x1*x1*x2 + x1*x1 + x1*x2*x2 - 2*x1 + 4;
h = x1*x1 + x2*x2 - 2;
g = 0.25*x1*x1 +0.75*x2*x2 -1;
%%% the gradient functions
gradf1 = diff(f,x1);
gradf2 = diff(f,x2);
% the hessian
hess = [diff(gradf1,x1), diff(gradf1,x2); diff(gradf2,x1), diff(gradf2,x2)];
% gradient of the constraints
gradh1 = diff(h,x1);
gradh2 = diff(h,x2);
gradg1 = diff(g,x1);
gradg2 = diff(g,x2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% graphical/symbolic solution for SLP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('***********************')
fprintf('\nSQP - Example 7.1')
fprintf('\n*********************\n')
for i = 1:it
%figure;
f1 = subs(f,{x1,x2},{xb(1),xb(2)});
g1 = subs(g,{x1,x2},{xb(1),xb(2)});
h1 = subs(h,{x1,x2},{xb(1),xb(2)});
%%% Print Information
fprintf('\n******************************')
fprintf('\nIteration : '),disp(i)
fprintf('******************************\n')
fprintf('Linearized about [x1, x2] : '),disp([xb(1) xb(2)])
fprintf('Objective function value f(x1,x2) : '),disp(f1);
fprintf('Equality constraint value value h(x1,x2) : '),disp(h1);
fprintf('Inequality constraint value value g(x1,x2) : '),disp(g1);
%fprintf('\nsolution for [delx1 delx2] : '),disp(sol')
% hold on
% calculate the value of the gradients
% f1 = subs(f,{x1,x2},{xb(1),xb(2)});
% g1 = subs(g,{x1,x2},{xb(1),xb(2)});
% h1 = subs(h,{x1,x2},{xb(1),xb(2)});
fprintf('\n-----------------------')
fprintf('\nQP - SUB PROBLEM')
fprintf('\n---------------------\n')
gf1 = double(subs(gradf1,{x1,x2},{xb(1),xb(2)}));
gf2 = double(subs(gradf2,{x1,x2},{xb(1),xb(2)}));
hess1 = double(subs(hess,{x1,x2},{xb(1),xb(2)}));
gh1 = double(subs(gradh1,{x1,x2},{xb(1),xb(2)}));
gh2 = double(subs(gradh2,{x1,x2},{xb(1),xb(2)}));
gg1 = double(subs(gradg1,{x1,x2},{xb(1),xb(2)}));
gg2 = double(subs(gradg2,{x1,x2},{xb(1),xb(2)}));
% the QP subproblem
syms dx1 dx2 % change in design
fquad = f1 + [gf1 gf2]*[dx1; dx2] + 0.5*[dx1 dx2]*hess1*[dx1 ;dx2];
hlin = h1 + [gh1 gh2]*[dx1; dx2];
glin = g1 + [gg1 gg2]*[dx1; dx2];
Fquadstr = strcat(num2str(f1),' + ',num2str(gf1), ...
'*','dx1',' + ',num2str(gf2),' * ','dx2', ...
' + 0.5*',num2str(hess1(1,1)),' * dx1^2', ...
' +',num2str(hess1(1,2)),' * dx1*dx2', ...
' + 0.5*',num2str(hess1(2,2)),' * dx2^2');
hlinstr = strcat(num2str(h1),' + ',num2str(gh1), ...
'*','dx1',' + ',num2str(gh2),' * ','dx2');
glinstr = strcat(num2str(g1),' + ',num2str(gg1), ...
'*','dx1',' + ',num2str(gg2),' * ','dx2');
fprintf('Quadratic Objective function f(x1,x2): \n'),disp(Fquadstr);
fprintf('Linearized equality h(x1,x2): '),disp(hlinstr);
fprintf('Linearized inequality g(x1,x2): '),disp(glinstr);
fprintf('\n')
% define Lagrangian for the QP problem
syms lamda beta
F = fquad + lamda*hlin + beta*glin;
fprintf('Case a: beta = 0\n');
Fnobeta = fquad + lamda*hlin;
%%% initialize best solution
dx1best = 0;
dx2best = 0;
Fbbest = 0;
%%%%%%%%%%%%%%%%%%%%%%%
%%% solve case (a) %%%
%%%%%%%%%%%%%%%%%%%%%%%
xcasea = solve(diff(Fnobeta,dx1),diff(Fnobeta,dx2),hlin);
sola = [double(xcasea.dx1) double(xcasea.dx2) double(xcasea.lamda)];
dx1a = double(xcasea.dx1);
dx2a = double(xcasea.dx2);
lamdaa = double(xcasea.lamda);
hlina = double(subs(hlin,{dx1,dx2},{dx1a,dx2a}));
glina = double(subs(glin,{dx1,dx2},{dx1a,dx2a}));
Fa = double(subs(Fnobeta,{dx1,dx2,lamda},{dx1a,dx2a,lamdaa}));
%%% results for case (a)
x1a = dx1a + xb(1);
x2a = dx2a + xb(2);
fv = double(subs(f,{x1,x2},{x1a,x2a}));
hv = double(subs(h,{x1,x2},{x1a,x2a}));
gv = double(subs(g,{x1,x2},{x1a,x2a}));
fprintf('Change in design vector: '),disp([dx1a dx2a]);
fprintf('The linearized quality constraint: '),disp(hlina);
fprintf('The linearized inequality constraint: '),disp(glina);
fprintf('New design vector: '),disp([x1a x2a]);
fprintf('The objective function: '),disp(fv);
fprintf('The equality constraint: '),disp(hv);
fprintf('The inequality constraint: '),disp(gv);
if (glina <= 0)
xb(1) = xb(1) + dx1a;
xb(2) = xb(2) + dx2a;
fbest = Fa;
dx1best = dx1a;
dx2best = dx2a;
end
%%%%%%%%%%%%%%%%%%%%%%%
%%% solve case (b) %%%
%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n Case b: g = 0\n');
xcaseb = solve(diff(F,dx1),diff(F,dx2),hlin,glin);
solb = [double(xcaseb.dx1) double(xcaseb.dx2) double(xcaseb.lamda) double(xcaseb.beta)];
dx1b = double(xcaseb.dx1);
dx2b = double(xcaseb.dx2);
betab = double(xcaseb.beta);
lamdab = double(xcaseb.lamda);
hlinb = double(subs(hlin,{dx1,dx2},{dx1b,dx2b}));
glinb = double(subs(glin,{dx1,dx2},{dx1b,dx2b}));
Fb = double(subs(F,{dx1,dx2,lamda,beta},{dx1b,dx2b,lamdab,betab}));
x1b = dx1b + xb(1);
x2b = dx2b + xb(2);
fv = double(subs(f,{x1,x2},{x1b,x2b}));
hv = double(subs(h,{x1,x2},{x1b,x2b}));
gv = double(subs(g,{x1,x2},{x1b,x2b}));
fprintf('Change in design vector: '),disp([dx1b dx2b]);
fprintf('The linearized quality constraint: '),disp(hlinb);
fprintf('The linearized inequality constraint: '),disp(glinb);
fprintf('New design vector: '),disp([x1b x2b]);
fprintf('The objective function: '),disp(fv);
fprintf('The equality constraint: '),disp(hv);
fprintf('The inequality constraint: '),disp(gv);
fprintf('Multiplier beta: '),disp(betab);
fprintf('Multiplier lamda: '),disp(lamdab);
if (betab > 0) & (Fb <= fbest)
xb(1) = x1b;
xb(2) = x2b;
dx1best = dx1b;
dx2best = dx2b;
end
%%% stopping criteria
if ([dx1best dx2best]*[dx1best dx2best]') <= 1.0e-08
fprintf('\n&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
fprintf('\nStopped: Design Not Changing')
fprintf('\n&&&&&&&&&&&&&&&&&&&&&&&&&&&&&\n\n')
break;
elseif i == it
fprintf('\n&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')
fprintf('\nStpped: Number of iterations at maximum')
fprintf('\n&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&\n\n')
break;
end
end
f1, g1, h1仍然是syms变量类型。- 在应用前使用函数
double()将它们更改为numeric类型 函数num2str() - 您已经对变量应用了
double()
gf1, gf2, gh1, gh2,gg1, gg2 and hess1
就在上面，不用碰它们
这是您应该替换的部分
Fquadstr = strcat(num2str(f1),' + ',num2str(gf1), ...
'*','dx1',' + ',num2str(gf2),' * ','dx2', ...
' + 0.5*',num2str(hess1(1,1)),' * dx1^2', ...
' +',num2str(hess1(1,2)),' * dx1*dx2', ...
' + 0.5*',num2str(hess1(2,2)),' * dx2^2');
hlinstr = strcat(num2str(h1),' + ',num2str(gh1), ...
'*','dx1',' + ',num2str(gh2),' * ','dx2');
glinstr = strcat(num2str(g1),' + ',num2str(gg1), ...
'*','dx1',' + ',num2str(gg2),' * ','dx2');
据此
% apply double() to f1
Fquadstr = strcat(num2str(double(f1)),' + ',num2str(gf1), ...
'*','dx1',' + ',num2str(gf2),' * ','dx2', ...
' + 0.5*',num2str(hess1(1,1)),' * dx1^2', ...
' +',num2str(hess1(1,2)),' * dx1*dx2', ...
' + 0.5*',num2str(hess1(2,2)),' * dx2^2');
% apply double() to h1
hlinstr = strcat(num2str(double(h1)),' + ',num2str(gh1), ...
'*','dx1',' + ',num2str(gh2),' * ','dx2');
% apply double() to g1
glinstr = strcat(num2str(double(g1)),' + ',num2str(gg1), ...
'*','dx1',' + ',num2str(gg2),' * ','dx2');