编写⽬标函数，定义放在⼀个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相似。
例：
matlab求导
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）
<后⾯是线性约束，⼤括号⾥⾯的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)
发布评论