matlab解决非线性规划问题(凸优化问题)

时间:2024-02-24 12:29:10

当目标函数含有非线性函数或者含有非线性约束的时候该规划问题变为非线性规划问题,非线性规划问题的最优解不一定在定义域的边界,可能在定义域内部,这点与线性规划不同;

例如:

 编写目标函数,定义放在一个m文件中;编写非线性约束条件函数矩阵,放在另一个m文件中;

function f = optf(x);
f = sum(x.^2)+8;
function [g, h] = limf(x);
g = [-x(1)^2+x(2)-x(3)^2
    x(1)+x(2)^2+x(3)^3-20]; %非线性不等式约束
h = [-x(1)-x(2)^2+2
     x(2)+2*x(3)^2-3];  %非线性等式约束
options = optimset(\'largescale\',\'off\');
[x y] = fmincon(\'optf\',rand(3,1),[],[],[],[],zeros(3,1),[],...
    \'limf\',options)

输出为:

最速下降法(求最小值):

代码如下:

function [f df] = detaf(x);
f = x(1)^2+25*x(2)^2;
df = [2*x(1)
      50*x(2)];
clc,clear;
x = [2;2];
[f0 g] = detaf(x);
while norm(g)>1e-6 %收敛条件为一阶导数趋近于0
    p = -g/norm(g);
    t = 1.0; %设置初始步长为1个单位
    f = detaf(x+t*p);
    while f>f0
        t = t/2;
        f = detaf(x+t*p);
    end %这一步很重要,为了保证最后收敛,保持f序列为一个单调递减的序列,否则很有可能在极值点两端来回震荡,最后无法收敛到最优值。
    x = x+t*p;
    [f0,g] = detaf(x);
end
x,f0

所得到的最优值为近似解。

第二种解法求极值是Newton法;

 

先写出ntfun.m和main.m,如下:

clc,clear;
x = [2;2];
[f0 g1 g2] = ntfun(x);
while norm(g1) > 1e-6 %收敛条件为一阶导趋于0,(如果二阶导为0那么p = 0,x会在鞍点处无法移动,所以newton法有一定缺陷,适合找驻点)我们一般用最速下降的方法来求最小值。
    p = -inv(g2)*g1;
    x = x + p;
    [f0,g1,g2] = ntfun(x);
end
x,f0
function [f, df d2f] = ntfun(x);
f = x(1)^4+25*x(2)^4+x(1)^2*x(2)^2;
df = [4*x(1)^3+2*x(1)*x(2)^2; 100*x(2)^3+2*x(1)^2*x(2)];
d2f = [12*x(1)^2+2*x(2)^2, 4*x(1)*x(2)
       4*x(1)*x(2), 300*x(2)^2+2*x(1)^2];
   

输出结果为近似值:

对于二次以上的函数,牛顿法并不能保证一定能找到最优点,只能找到局部最优点或者鞍点,

我们可以仿照梯度下降方法改进牛顿法。

clc,clear
x=[2;2]; [f0,g1,g2]=ntfun(x); 
while norm(g1)>1e-5
p=-inv(g2)*g1;
p=p/norm(p); 
t=1.0;
f=ntfun(x+t*p); 
while f>f0
t=t/2;
f=ntfun(x+t*p);
end
x=x+t*p; 
[f0,g1,g2]=ntfun(x); 
end
x,f0

输出为:

 变尺度法(变步长):

常用Hesse矩阵来逼近二阶导数矩阵;

 

 

这种方法相对于牛顿法和梯度下降法收敛速度更快,但是计算量相对较大较复杂,并且只适用于可以求导或者近似求导的问题。

直接法(Powell方法):

 

 Matlab 无约束求极值问题:

matlab有两个函数,一个是fminunc(fun,x0,options,p1,p2,..) 当fun只有一个返回值时候,返回函数是f(x);如果有三个返回值,则第二个为梯度函数向量,第三个为二阶导数矩阵。options 为优化参数,可以使用缺省参数,p1,p2是传递给fun的一些参数。

例:

options = optimset(\'GradObj\',\'on\');
[x,y] = fminunc(\'fun\',rand(1,2),options)
function [f df] = fun(x);
f = 100*(x(2)-x(1)^2)^2+(1-x(1))^2;
df = [-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)];

如果要加Hessian矩阵(函数的二阶导矩阵),则main代码如下,fun中应该第三项返回d2f的函数矩阵:

options = optimset(\'GradObj\',\'on\',\'Hessian\',\'on\');
[x,y] = fminunc(\'fun\',rand(1,2),options)

 fminsearch函数可以求得初值附近的极小点和极小值,用法与fminunc相似。

 

例:

H=[4,-4;-4,8];
f=[-6;-3];
a=[1,1;4,1];
b=[3;9]; 
[x,value]=quadprog(H,f,a,b,[],[],zeros(2,1))

输出:

 

惩罚函数法:

 

 

例:

求解程序如下:

function g = fun(x);
M = 1e5; %M充分大
f = x(1)^2+x(2)^2+8; %目标函数
g = f-M*min(x(1),0)-M*min(x(2),0)-M*min(x(1)^2-x(2),0)+...
    M*abs(-x(1)-x(2)^2+2); %构造的新目标函数
[x y] = fminunc(\'fun\',rand(1,2))

输出:

Matlab求约束极值问题相关的函数:

 

单变量非线性函数在区间上的极小值: fminbnd函数

 fseminf函数:半无穷约束(sem + inf)

 s.t.后面是线性约束,大括号里面的C(x)为不等式非线性约束,Ceq(x)为等式非线性约束,PHI(x,w)为附加变量w对x的约束,

注意F(x)中可是没有变量w的。PHI(x,w)<=0为半无穷约束,所以这个函数也叫做fseminf.

NTHETA是半无穷约束PHI(x,w)的个数,下题中该值为2,k1&k2;

 

function f = fun(x,s);
f = sum((x-0.5).^2);
function [c ceq k1 k2 s] = fun2(x,s) %定义seminfuncon s表示推荐的取样步长
c = [];
ceq = [];
if isnan(s(1,1))
    s = [0.2 0;0.2 0];
end
%取样值
w1 = 1:s(1,1):100;
w2 = 1:s(2,1):100;
%半无穷约束
k1=sin(w1*x(1)).*cos(w1*x(2))-1/1000*(w1-50).^2-sin(w1*x(3))-x(3)-1; 
k2=sin(w2*x(2)).*cos(w2*x(1))-1/1000*(w2-50).^2-sin(w2*x(3))-x(3)-1;
%画出半无穷约束的图形
plot(w1,k1,\'-\',w2,k2,\'+\');
[x y] = fseminf(@fun,rand(3,1),2,@fun2)

  由图可见对于任意的w1,w2,对应的k1,k2均小于0,满足了这两个半无穷约束。

程序通过依据步长s遍历每一组(w1,w2),将其转为许多个非线性不等式约束,即将fseminfuncon转为许多个C(x),然后综合在一起,

再求得在这接近无数组的非线性不等式约束下的极值点。

fminimax函数(fmin imax)求上界函数的最小值点的函数,max Fi(x) 表示一族函数集

中最大的那个函数

function f = fun(x);
f = [2*x(1)^2+x(2)^2-48*x(1)-40*x(2)+304
     -x(1)^2-3*x(2)^2
     x(1)+3*x(2)-18
     -x(1)-x(2)
     x(1)+x(2)-8];
[x y] = fminimax(@fun,rand(2,1))

输出:

 在选择用梯度优化选项时候,除了f需要给出df之外,非线性约束函数也要给出对应的梯度,例如dc和dceq,

按行求导,第一行对x1进行求导,第二行对x2进行求导。

非线性规划这一章节结束。