當目標函數含有非線性函數或者含有非線性約束的時候該規划問題變為非線性規划問題,非線性規划問題的最優解不一定在定義域的邊界,可能在定義域內部,這點與線性規划不同;
例如:
編寫目標函數,定義放在一個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進行求導。
非線性規划這一章節結束。