數學規划模型——非線性規划問題

本文轉載自查看原文 2020-02-27 19:41 1468 數學建模

title: 數學規划模型——非線性規划問題
date: 2020-02-26 20:10:37
categories: 數學建模
tags: [MATLAB, 數學規划模型]

Matlab 中⾮線性規划的標准型

\[min \ f(x) \\ \ \\ s.t. \ \begin{cases} & \text{} AX<=b,Aeg*x=beg \quad 線性約束\\ & \text{} c(x)<=0,ceg(x)=0\quad 非線性約束 \\ & \text{}lb<=x<=ub \quad 上下界約束(也可以當成不等式約束) \end{cases} \]

練習題

Matlab求解⾮線性規划的函數

[x,fval] = fmincon(@fun,x0,A,b,Aeq,beq,lb,ub,@nonlfun,option)

⾮線性規划中對於初始值 x0 的選取⾮常重要 , 因為⾮線性規划的算法求解出來的是⼀個局部優解。 (線性規划不存在這個問題 )
如果要求全局優解! 有兩種思路 :
- 給定不同的初始值 , 在⾥⾯找到⼀個優解 ;
- 先⽤蒙特卡羅模擬, 得到⼀個蒙特卡羅解 , 然后將這個解作為初始值來求優解
"option" 選項可以給定求解的算法, ⼀共有四種 :
- interior.point (內點法)
- sqp(序列⼆次規划法)
- active.set (有效集法)
- trust-region.reflective (信賴域反射算法) 。
不同的算法有其各⾃的優缺點和適⽤情況 , 我們可以改變求解的算法來看求解的結果是否變好了。如何改變求解的算法請參考代碼演示。 (數模⽐賽中可以都嘗試下, 這可以體現出穩健性, 也是你的優點)
"@fun" 表示⽬標函數
我們要編寫⼀個獨⽴的 m ⽂件來存⽬標函數 :
```
function f = fun ( x) 
    f=...
end
```
- 注 fun 可以任意取名 , 只要符合 Matlab 命名規范 , 保存的m ⽂件也是這個名。
- f 也可以任意取名, 但返回的 f 和函數內部的f得完全⼀致;
- 這⾥的x實際上是表示決策變量的向量 ,其行列方向取決於初始值。
- 調⽤函數 : fmincon(@fun,...) 求解。
" @nonlfun " 表示⾮線性部分的約束 , 我們同樣得編寫⼀個獨⽴的 m ⽂件儲存非線性約束條件 :
- nonlfun 同樣可任意取名 , 別和上⾯的fun相同即可, 保存的m文件也得是這個名字。
- c,ceq中可能有多個約束因此我們寫成列向量的形式 ;
- 若不存在⾮線性不等式約束 , 則可以令 C = []
- 調⽤函數 fmincon(.......,@nonlfun.options) 求解。
注意要把下標改寫為擴號 , 例如\(f=x_{1}^2+3x_{2}\)名寫成 Matlab 能識別的就應該為 : \(f=x(1)^2+3*x(2)\)
若不存在某種約束 , 則可⽤ " [] "替代 , 若后⾯全為 "[ ] " 且不指定option(使⽤默認的求解⽅法) , 則"[]"也可以省略掉。

練習題代碼

code.m

% 非線性規划的函數
% [x,fval] = fmincon(@fun,x0,A,b,Aeq,beq,lb,ub,@nonlfun,option)
% x0表示給定的初始值（用行向量或者列向量表示），必須得寫
% A b表示線性不等式約束
% Aeq beq 表示線性等式約束
% lb ub 表示上下界約束
% @fun表示目標函數
% @nonlfun表示非線性約束的函數
% option 表示求解非線性規划使用的方法
clear;clc
format long g   %可以將Matlab的計算結果顯示為一般的長數字格式（默認會保留四位小數，或使用科學計數法）

%% 例題1的求解
% max f(x) = x1^2 +x2^2 -x1*x2 -2x1 -5x2
% s.t. -(x1-1)^2 +x2 >= 0 ;  2x1-3x2+6 >= 0
x0 = [0 0];  %任意給定一個初始值 
A = [-2 3]; b = 6;
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1)  % 注意 fun1.m文件和nonfun1.m文件都必須在當前文件夾目錄下
fval = -fval
% 一個值得討論的地方，能不能把線性不等式約束Ax <= b也寫到nonlfun1函數中？
% 先把nonlfun1中的c改為下面這樣：
% c = [(x(1)-1)^2-x(2); 
%        -2*x(1)+3*x(2)-6];
%  [x,fval] = fmincon(@fun1,x0,[],[],[],[],[],[],@nonlfun1)
% 結果也是可以計算出來的，但並不推薦這樣做~

%% 使用其他算法對例題1求解
% edit fmincon  % 查看fmincon的“源代碼”
% Matlab2017a默認使用的算法是'interior-point' 內點法
% 使用interior point算法 （內點法）
option = optimoptions('fmincon','Algorithm','interior-point')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)  
fval = -fval
% 使用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  
% 使用active set算法 （有效集法）
option = optimoptions('fmincon','Algorithm','active-set')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)
fval = -fval  
% 使用trust region reflective (信賴域反射算法)
option = optimoptions('fmincon','Algorithm','trust-region-reflective')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)  
fval = -fval  
% this algorithm does not solve problems with the constraints you have specified. 
% 這說明這個算法不適用我們這個約束條件，所以以后遇到了不能求解的情況，記得更換其他算法試試！！！

%% 選取初始值得到的結果可能會不滿足限定條件，出現了一個Bug 因此選擇的初始值很重要
x0 = [40.8, 10.8];
option = optimoptions('fmincon','Algorithm','interior-point')
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option)  
fval = -fval  
% https://cn.mathworks.com/help/optim/ug/fmincon.html

%% 生成不同的隨機初始值來優化代碼，有一定幾率會觸發上面那個Bug，因此不推薦
n = 10;  % 重復n次
Fval = +inf; X = [0,0];  %初始化最優的結果
A = [-2 3]; b = 6;
for i = 1:n
    x0 = [rand()*10 , rand()*10];  %用隨機數生成一個初始值（隨機數的范圍自己根據題目條件設置） 
    [x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option); % 注意 fun1.m文件和nonfun1.m文件都必須在當前文件夾目錄下
    if fval < Fval  % 如果找到了更小的值，那么就代替最優的結果
        Fval = fval;
        X = x;
    end
end
Fval = -Fval
X

%% 使用蒙特卡羅的方法來找初始值(推薦）
clc,clear;
n=10000000; %生成的隨機數組數
x1=unifrnd(-100,100,n,1);  % 生成在[-100,100]之間均勻分布的隨機數組成的n行1列的向量構成x1
x2=unifrnd(-100,100,n,1);  % 生成在[-100,100]之間均勻分布的隨機數組成的n行1列的向量構成x2
fmin=+inf; % 初始化函數f的最小值為正無窮（后續只要找到一個比它小的我們就對其更新）
for i=1:n
    x = [x1(i), x2(i)];  %構造x向量, 這里千萬別寫成了：x =[x1, x2]
    if ((x(1)-1)^2-x(2)<=0)  & (-2*x1(i)+3*x2(i)-6 <= 0)     % 判斷是否滿足條件
        result = -x(1)^2-x(2)^2 +x(1)*x(2)+2*x(1)+5*x(2) ;  % 如果滿足條件就計算函數值
        if  result  < fmin  % 如果這個函數值小於我們之前計算出來的最小值
            fmin = result;  % 那么就更新這個函數值為新的最小值
            x0 = x;  % 並且將此時的x1 x2更新為初始值
        end
    end
end
disp('蒙特卡羅選取的初始值為：'); disp(x0)
A = [-2 3]; b = 6;
[x,fval] = fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1)
fval = -fval  



%% 例題二的求解
x0 = [1 1 1];  %任意給定一個初始值 
lb = [0 0 0];  % 決策變量的下界
[x,fval] = fmincon(@fun2,x0,[],[],[],[],lb,[],@nonlfun2)  % 注意 fun2.m文件和nonfun2.m文件都必須在當前文件夾目錄下
% x =
%          0.552167405729277          1.20325915507969         0.947824046150443
% fval =
%           10.6510918606939



%% 使用蒙特卡羅的方法來找初始值(推薦）
clc,clear;
n=1000000; %生成的隨機數組數
x1= unifrnd(0,2,n,1);   % 生成在[0,2]之間均勻分布的隨機數組成的n行1列的向量構成x1
x2 = sqrt(2-x1);  % 根據非線性等式約束用x1計算出x2
x3 = sqrt((3-x2)/2); % 根據非線性等式約束用x2計算出x3
fmin=+inf; % 初始化函數f的最小值為正無窮（后續只要找到一個比它小的我們就對其更新）
for i=1:n
    x = [x1(i), x2(i), x3(i)];  %構造x向量, 這里千萬別寫成了：x =[x1, x2, x3]
    if (-x(1)^2+x(2)-x(3)^2<=0) & (x(1)+x(2)^2+x(3)^2-20<=0)   % 判斷是否滿足條件
        result =sum(x.*x) + 8 ;  % 如果滿足條件就計算函數值
        if  result  < fmin  % 如果這個函數值小於我們之前計算出來的最小值
            fmin = result;  % 那么就更新這個函數值為新的最小值
            x0 = x;  % 並且將此時的x1 x2 x3更新為初始值
        end
    end
end
disp('蒙特卡羅選取的初始值為：'); disp(x0)
lb = [0 0 0];  % 決策變量的下界
[x,fval] = fmincon(@fun2,x0,[],[],[],[],lb,[],@nonlfun2)  % 注意 fun2.m文件和nonfun2.m文件都必須在當前文件夾目錄下

%% 例題三的求解(蒙特卡羅模擬那一講的例題)
clear;clc
% 蒙特卡羅模擬得到的最大值為3445.6014
% 最大值處x1 x2 x3的取值為：
%           22.5823101903968          12.5823101903968          12.1265223966757
A = [1 -2 -2;  1 2 2];  b = [0 72];
x0 = [ 22.58   12.58  12.13];
Aeq = [1 -1 0]; beq = 10;
lb = [-inf 10 -inf];  ub = [inf 20 inf];  
[x,fval] = fmincon(@fun3,x0,A,b,Aeq,beq,lb,ub,[])  % 注意沒有非線性約束，所以這里可以用[]替代，或者干脆不寫
fval = -fval

fun1.m

function f = fun1(x)
    % 注意：這里的f實際上就是目標函數，函數的返回值也是f
    % 輸入值x實際上就是決策變量，由x1和x2組成的向量
    % fun1是函數名稱，到時候會被fmincon函數調用, 可以任意取名
    % 保存的m文件和函數名稱得一致，也要為fun1.m
%      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

fun2.m

function f = fun2(x)
    %     f = x(1)^2+x(2)^2 +x(3)^2+8 ; 
    f = sum(x.*x) + 8;  % 可別忘了x實際上是一個向量，我們可以使用矩陣的運算符號對其計算
end

fun3.m

function f = fun3(x)
    f = -prod(x);  % 可別忘了x實際上是一個向量（prod表示連乘符號，用法和sum類似）
end

nonlfun1.m

function [c,ceq] = nonlfun1(x)
    % 注意：這里的c實際上就是非線性不等式約束，ceq實際上就是非線性等式約束
    % 輸入值x實際上就是決策變量，由x1和x2組成的一個向量
    % 返回值有兩個，一個是非線性不等式約束c，一個是非線性等式約束ceq
    % nonlfun1是函數名稱，到時候會被fmincon函數調用, 可以任意取名，但不能和目標函數fun1重名
    % 保存的m文件和函數名稱得一致，也要為nonlfun1.m
%     -(x1-1)^2 +x2 >= 0 
   c = [(x(1)-1)^2-x(2)];   % 千萬別寫成了: (x1-1)^2 -x2
   ceq = [];  % 不存在非線性等式約束，所以用[]表示
end

nonlfun2.m

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

典型例題

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 數學規划模型——線性規划問題 Python數學建模系列（三）：規划問題之非線性規划 Matlab 非線性規划問題模型代碼非線性規划非線性規划（1）線性規划和非線性規划 [數學建模(五)]線性規划，二次規划和非線性規划的MATLAB求解 Matlab非線性規划 MATLAB規划問題——線性規划和非線性規划(from CSDN echo_hello) MATLAB 非線性規划