matlab求解非線性規劃

卢宇博發表於2024-08-11
Matlab

目錄
  • 前言
  • 一、非線性規劃的標準型
  • 二、fmincon函式
    • 1.目標函式--function f = fun(x)
    • 2.非線性約束函式--[c,ceq] = nonlfun(x)
    • 3.設定求解方法--option
  • 三、matlab求解非線性規劃的例項與可能遇到的問題
    • 1.初值問題
    • 2.演算法問題
      • （1）內點法求解
      • （2）SQP演算法求解
      • （3）active set演算法求解
      • （4）信賴域反射演算法求解

前言

描述目標函式或約束條件的數學表示式中,至少一個是非線性函式，這樣的最佳化問題通常稱為非線性規劃。一般說來,解決非線性規問題要比解決線性規劃問題困難得多。線性規劃有適用於一般情況的單純形法,但是於非線性規劃問題,目前還沒有一種適用於一般情況的求解方法,現有各種方法都有各特定的適用範圍。本章會介紹如何用matlab內建的不同演算法解決非線性規劃問題。

一、非線性規劃的標準型

非線性規劃的標準型相比線性規劃多了一個非線性不等式約束與非線性等式約束，其與線性約束的區別在於非線性約束需要把等號的右邊全部化成0

二、fmincon函式

[x,fval] = fmincon(@fun,x0,A,b,Aeq,beq,lb,ub,@nonlfun,option)
x0表示給定的初始值（用行向量或者列向量表示），必須得寫，而且對求解十分重要
@fun表示目標函式，不再像線性規劃的函式僅用係數向量這麼簡單來表示，需要在同一檔案目錄下建立一個函式檔案
@nonlfun表示非線性約束的函式
option 表示求解非線性規劃使用的方法

1.目標函式--function f = fun(x)

此引數是一個儲存在同一資料夾下的函式檔案，其中x是一個向量，寫的時候不能直接寫x1,x2；而是要用x(1)，x(2)來表示向量中的元素

function f = fun(x)
%      max  f(x) = x1^2 +x2^2 -x1*x2 -2x1 -5x2
    f = -x(1)^2-x(2)^2 +x(1)*x(2)+2*x(1)+5*x(2) ; 
end

2.非線性約束函式--[c,ceq] = nonlfun(x)

這裡的c是一個非線性不等式約束，ceq是非線性等式約束，他們的寫法也是把x當作一個向量，用x(1)來引用x裡面的元素。
同時，需要把約束的右側全部化為0，填入約束的左側

function [c,ceq] = nonlfun2(x)
    % 非線性不等式約束
    c = [-x(1)^2+x(2)-x(3)^2;   % 一定要注意寫法的規範，再次強調這裡的x是一個向量！不能把x(1)寫成x1
            x(1)+x(2)^2+x(3)^2-20];
    % 非線性等式約束
    ceq = [-x(1)-x(2)^2+2;
                x(2)+2*x(3)^2-3]; 
end

3.設定求解方法--option

matlab為求解非線性約束提供了幾種方法，可以透過命令來呼叫，在比賽中，我們可以在論文中說明我們使用了不同方法，並描述哪一種方法更好
（1）內點法-- matlab預設的求解方法
option = optimoptions('fmincon','Algorithm','interior-point')
（2）序列二次規劃法
option = optimoptions('fmincon','Algorithm','sqp')
（3）有效集法
option = optimoptions('fmincon','Algorithm','active-set')
（4）信賴域反射演算法（對約束有要求）
option = optimoptions('fmincon','Algorithm','trust-region-reflective')

三、matlab求解非線性規劃的例項與可能遇到的問題

1.初值問題

x0 = [0 0]; x1 = [40.8, 10.8] %任意給定一個初始值 
A = [-2 3]; b = 6;
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1)  % 注意 fun1.m檔案和nonlfun1.m檔案都必須在當前資料夾目錄下
fval = -fval

當選定x0作為初值時，結果：x=[1.00000001604465,6.66655111798166e-08]
fval=-1.00000039999306
當選定x1作為初值時，結果：x=[1.5123164061584e+22,-1.85553223470542e+52]
fval=3.4429998740309e+104
這說明初值對求解的影響很大，不同的初值給出的結果天差地別，而且讓x1作為初值時，得到的結果應該也不滿足我們的約束。（或許我們的初值就不滿足約束）
因此，找到一個合適的初值對解決我們的問題時非常重要的
使用蒙特卡羅的方法來找初始值(推薦），在蒙特卡羅的方法的迴圈過程中，不必再加入判斷等式約束的條件，因為我們可以透過等式約束的條件來生成隨機數

2.演算法問題

（1）內點法求解

option = optimoptions('fmincon','Algorithm','interior-point')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)  
fval = -fval

x =[1.00000001604465,6.66655111798166e-08]
fval =-1.00000039999306

（2）SQP演算法求解

option = optimoptions('fmincon','Algorithm','sqp')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)  
fval = -fval   %得到-4.358，遠遠大於內點法得到的-1,猜想是初始值的影響
% 改變初始值試試
x0 = [1 1];  %任意給定一個初始值 
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)  % 最小值為-1，和內點法相同（這說明內點法的適應性要好）
fval = -fval

（3）active set演算法求解

option = optimoptions('fmincon','Algorithm','active-set')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)
fval = -fval

結果：
x =[-0.333333399760948,1.7777777334927]
fval =-4.35802434691854

（4）信賴域反射演算法求解

option = optimoptions('fmincon','Algorithm','trust-region-reflective')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)  
fval = -fval


這說明這個演算法不適用我們這個約束條件，所以以後遇到了不能求解的情況，記得更換其他演算法試試！！！

相關文章