matlab解决非线性规划问题（凸优化问题）

matlab解决⾮线性规划问题（凸优化问题）
当⽬标函数含有⾮线性函数或者含有⾮线性约束的时候该规划问题变为⾮线性规划问题，⾮线性规划问题的最优解不⼀定在定义域的边界，可能在定义域内部，这点与线性规划不同；例如：
编写⽬标函数，定义放在⼀个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进⾏求导。

⾮线性规划这⼀章节结束。