劃實(shí)戰(zhàn)：從fmincon算法選擇到全局優(yōu)化策略)
1. 項(xiàng)目概述從線性到非線性的思維躍遷在數(shù)學(xué)建模的實(shí)戰(zhàn)中我們遇到的絕大多數(shù)問題其目標(biāo)函數(shù)或約束條件都不是簡(jiǎn)單的線性關(guān)系。比如你想優(yōu)化一個(gè)工廠的生產(chǎn)計(jì)劃成本可能隨著產(chǎn)量呈指數(shù)增長(zhǎng)目標(biāo)函數(shù)非線性或者你設(shè)計(jì)一個(gè)機(jī)械結(jié)構(gòu)其應(yīng)力必須小于材料的非線性屈服強(qiáng)度約束條件非線性。這時(shí)線性規(guī)劃那套漂亮的單純形法就完全失效了。非線性規(guī)劃正是為了解決這類“彎彎繞繞”的優(yōu)化問題而生的核心數(shù)學(xué)工具。它不像線性規(guī)劃那樣有“標(biāo)準(zhǔn)答案”式的通用解法更像是一個(gè)工具箱里面裝著各種針對(duì)不同問題特性的“專用扳手”。我接觸過很多剛開始做建模的同學(xué)一看到“非線性”三個(gè)字就頭疼覺得深不可測(cè)。其實(shí)不然它的核心思想非常直觀在復(fù)雜的地形目標(biāo)函數(shù)曲面上找到那個(gè)最低點(diǎn)最小值或最高點(diǎn)最大值同時(shí)不能跑到禁區(qū)約束條件里去。這次我們就來徹底拆解這個(gè)工具箱不僅告訴你每個(gè)工具算法怎么用更重點(diǎn)講清楚什么時(shí)候該用哪個(gè)以及用的時(shí)候最容易在哪兒翻車。我們會(huì)以最常用的MATLAB環(huán)境為例手把手帶你從理論走到代碼實(shí)現(xiàn)讓你下次遇到非線性問題時(shí)能胸有成竹地選出最合適的那把“扳手”。2. 非線性規(guī)劃的核心思想與問題分類在動(dòng)手寫代碼之前我們必須先搞清楚面對(duì)的是什么“型號(hào)”的問題。非線性規(guī)劃問題通常可以寫成如下標(biāo)準(zhǔn)形式最小化問題Minimize: f(x) Subject to: g_i(x) ≤ 0, i 1, ..., m (不等式約束) h_j(x) 0, j 1, ..., p (等式約束) x ∈ R^n (決策變量)這里f(x), g_i(x), h_j(x) 中至少有一個(gè)是非線性函數(shù)。根據(jù)這些函數(shù)的特性我們可以把問題分門別類這直接決定了我們?cè)撨x用哪種算法。2.1 凸與非凸決定問題難度的分水嶺這是非線性規(guī)劃中最關(guān)鍵的分類沒有之一。凸規(guī)劃如果目標(biāo)函數(shù) f(x) 是凸函數(shù)并且不等式約束函數(shù) g_i(x) 是凸函數(shù)等式約束 h_j(x) 是線性函數(shù)那么這個(gè)問題就是凸規(guī)劃。凸規(guī)劃的任何局部最優(yōu)解必定是全局最優(yōu)解。這是它最大的優(yōu)點(diǎn)意味著算法只要找到一個(gè)“坑底”那就是整個(gè)區(qū)域的最低點(diǎn)。求解凸規(guī)劃相對(duì)“友好”。非凸規(guī)劃不滿足上述凸性條件的規(guī)劃問題。它的“地形圖”可能像連綿的群山有無數(shù)個(gè)山谷局部最優(yōu)點(diǎn)算法很容易陷在某個(gè)小山谷里而找不到最深的那一個(gè)全局最優(yōu)點(diǎn)。求解非凸規(guī)劃是NP-Hard問題通常只能尋找“較好的”局部最優(yōu)解或者采用一些隨機(jī)策略如模擬退火、遺傳算法來嘗試尋找全局最優(yōu)。實(shí)操心得在實(shí)際建模中我們首先應(yīng)該嘗試判斷問題是否具有凸性。一個(gè)簡(jiǎn)單的技巧如果目標(biāo)函數(shù)是二次型且Hessian矩陣半正定或者約束是線性的那么它很可能是凸的。對(duì)于復(fù)雜函數(shù)判斷凸性需要利用二階條件Hessian矩陣處處半正定這在實(shí)踐中往往很困難。因此一個(gè)務(wù)實(shí)的做法是默認(rèn)問題是非凸的然后選擇能處理非凸問題的穩(wěn)健算法同時(shí)嘗試從多個(gè)不同的初始點(diǎn)出發(fā)求解以降低陷入糟糕局部最優(yōu)的風(fēng)險(xiǎn)。2.2 無約束與有約束解決問題的基本框架無約束非線性優(yōu)化問題中沒有任何 g_i(x) 和 h_j(x) 的限制。這類問題的經(jīng)典算法構(gòu)成了非線性優(yōu)化的基石例如梯度下降法沿著目標(biāo)函數(shù)負(fù)梯度方向迭代簡(jiǎn)單但收斂慢。牛頓法利用目標(biāo)函數(shù)的二階導(dǎo)數(shù)Hessian矩陣信息收斂速度快但需要計(jì)算Hessian矩陣及其逆計(jì)算量大。擬牛頓法如BFGS, DFP通過構(gòu)造一個(gè)近似矩陣來模擬Hessian矩陣的逆既保持了較快的收斂速度又避免了直接計(jì)算Hessian矩陣是實(shí)踐中無約束優(yōu)化的首選。有約束非線性優(yōu)化這是我們討論的重點(diǎn)也是fmincon等求解器主要應(yīng)對(duì)的場(chǎng)景。核心思路是將有約束問題轉(zhuǎn)化為一系列無約束或更簡(jiǎn)單的約束問題來求解。2.3 二次規(guī)劃非線性中的“線性”特例二次規(guī)劃是指目標(biāo)函數(shù)是二次函數(shù)約束條件是線性函數(shù)的一類特殊非線性規(guī)劃。它的標(biāo)準(zhǔn)形式為 Minimize: (1/2) * x^T * H * x c^T * x Subject to: A * x ≤ b, Aeq * x beq, lb ≤ x ≤ ub雖然目標(biāo)函數(shù)是非線性的二次但由于其結(jié)構(gòu)特殊存在非常高效和可靠的專用算法如有效集法、內(nèi)點(diǎn)法。在MATLAB中可以使用quadprog函數(shù)專門求解QP問題。很多復(fù)雜的非線性問題在局部可以用二次函數(shù)來近似因此QP求解器也是許多高級(jí)非線性算法如序列二次規(guī)劃SQP的核心子步驟。3. MATLAB實(shí)戰(zhàn)核心fmincon求解器深度解析MATLAB的fmincon是求解中小規(guī)模有約束非線性規(guī)劃問題的“瑞士軍刀”。它的強(qiáng)大之處在于內(nèi)部集成了多種算法可以自動(dòng)或手動(dòng)選擇以適應(yīng)不同問題。但要用好它必須理解其每一個(gè)參數(shù)背后的意義。3.1fmincon的基本調(diào)用與參數(shù)精講一個(gè)最基礎(chǔ)的調(diào)用格式如下[x, fval, exitflag, output] fmincon(fun, x0, A, b, Aeq, beq, lb, ub, nonlcon, options)我們來逐一拆解這些輸入輸出參數(shù)fun目標(biāo)函數(shù)句柄。例如(x) x(1)^2 x(2)^2。這里最容易出錯(cuò)的地方是函數(shù)定義必須能接受向量輸入x并返回標(biāo)量值。如果目標(biāo)函數(shù)計(jì)算量很大可以考慮在函數(shù)內(nèi)部進(jìn)行向量化操作或使用全局變量/嵌套函數(shù)傳遞額外參數(shù)。x0初始猜測(cè)值。這是影響求解結(jié)果最關(guān)鍵的因素之一尤其對(duì)于非凸問題。糟糕的初始點(diǎn)可能導(dǎo)致算法收斂到很差的局部最優(yōu)甚至失敗。一個(gè)好的策略是根據(jù)物理意義或經(jīng)驗(yàn)給出初始值或者進(jìn)行簡(jiǎn)單的網(wǎng)格搜索從多個(gè)初始點(diǎn)中選取最好的結(jié)果。A, b, Aeq, beq, lb, ub線性約束和邊界約束。這是定義約束最高效的方式應(yīng)優(yōu)先使用。nonlcon非線性約束函數(shù)句柄。該函數(shù)需要返回兩個(gè)輸出[c, ceq]其中c(x) 0表示非線性不等式約束ceq(x) 0表示非線性等式約束。即使只有一種約束也必須同時(shí)返回兩個(gè)輸出將不存在的那個(gè)設(shè)為空數(shù)組[]。3.2 算法選擇options的設(shè)置藝術(shù)通過optimoptions(‘fmincon’)來設(shè)置選項(xiàng)。算法選擇 (Algorithm) 是核心‘interior-point’(內(nèi)點(diǎn)法)默認(rèn)且最通用的算法。特別適合大規(guī)模問題能高效處理邊界約束和稀疏性。它通過在可行域內(nèi)部構(gòu)造一條路徑逼近最優(yōu)解。對(duì)于大多數(shù)問題首選這個(gè)算法?！畇qp’(序列二次規(guī)劃)另一種強(qiáng)大的通用算法。它通過在每一步迭代中求解一個(gè)二次規(guī)劃子問題來尋找搜索方向。對(duì)于中小規(guī)模問題尤其是約束較多的問題表現(xiàn)可能比內(nèi)點(diǎn)法更好?！產(chǎn)ctive-set’(有效集法)一種較老的算法適用于中小規(guī)模問題。它能精確識(shí)別在最優(yōu)解處起作用的約束active constraints。對(duì)于需要知道哪些約束是“緊”的問題這個(gè)算法有優(yōu)勢(shì)?！畉rust-region-reflective’(信賴域反射法)這個(gè)算法要求目標(biāo)函數(shù)是非線性的標(biāo)量函數(shù)且只能處理邊界約束或線性等式約束不能直接處理非線性約束或線性不等式約束但可以通過轉(zhuǎn)換。它的優(yōu)勢(shì)在于能利用目標(biāo)函數(shù)的梯度信息對(duì)于特定類型的問題非常高效。避坑指南如果你的問題包含非線性約束那么算法只能從‘interior-point’和‘sqp’中選擇。初次求解一個(gè)未知問題時(shí)建議先使用默認(rèn)的‘interior-point’算法。如果收斂速度慢或不穩(wěn)定再嘗試‘sqp’。務(wù)必在選項(xiàng)中打開梯度檢查options optimoptions(‘fmincon’, ‘CheckGradients’, true);這能幫你發(fā)現(xiàn)自定義梯度函數(shù)中的錯(cuò)誤避免因梯度不準(zhǔn)導(dǎo)致算法失敗。3.3 輸出結(jié)果解讀與診斷求解完成后不能只看最優(yōu)解x和最優(yōu)值fvalexitflag和output包含了至關(guān)重要的診斷信息。exitflag(退出標(biāo)志) 0算法收斂到局部最優(yōu)解。這是成功標(biāo)志。 0迭代次數(shù)或函數(shù)計(jì)算次數(shù)超過了MaxIterations或MaxFunctionEvaluations選項(xiàng)設(shè)置的最大值。此時(shí)得到的解可能不是最優(yōu)的需要增加迭代上限或檢查問題 formulation。 0求解失敗。常見原因包括目標(biāo)函數(shù)或約束函數(shù)在迭代點(diǎn)處返回了NaN或Inf問題可能無界初始點(diǎn)不可行對(duì)于嚴(yán)格要求可行性的算法。需要根據(jù)output.message中的信息進(jìn)行排查。output結(jié)構(gòu)體包含迭代次數(shù) (iterations)、函數(shù)計(jì)算次數(shù) (funcCount)、一階最優(yōu)性條件 (firstorderopt)、算法類型等。firstorderopt衡量了當(dāng)前解滿足一階最優(yōu)性條件KKT條件的程度這個(gè)值越小說明解越“優(yōu)”。通常小于1e-6可以認(rèn)為是很好的收斂。4. 進(jìn)階技術(shù)與實(shí)戰(zhàn)策略掌握了fmincon的基本用法我們來看看如何解決更復(fù)雜的情況以及如何提升求解的效率和穩(wěn)定性。4.1 處理復(fù)雜非線性約束與可行性當(dāng)非線性約束非常復(fù)雜時(shí)算法可能很難找到一個(gè)可行的初始點(diǎn)或者在迭代中保持可行性。這時(shí)可以嘗試使用罰函數(shù)法這是將約束問題轉(zhuǎn)化為無約束問題的經(jīng)典思路?；舅枷胧菍⒓s束違反的程度作為一個(gè)“懲罰項(xiàng)”加到目標(biāo)函數(shù)中。例如對(duì)于約束g(x) 0可以構(gòu)造罰函數(shù)P(x) f(x) μ * max(0, g(x))^2其中μ是一個(gè)很大的正數(shù)罰因子。然后使用無約束優(yōu)化方法如fminunc求解P(x)。隨著μ增大解會(huì)越來越逼近原約束問題的解。缺點(diǎn)是罰因子需要精心選擇太大可能導(dǎo)致數(shù)值問題太小則約束得不到滿足。fmincon的可行性模式對(duì)于某些算法如sqp可以設(shè)置options.ConstraintTolerance來放寬對(duì)約束的嚴(yán)格滿足要求讓算法先找到一個(gè)“差不多”可行的點(diǎn)再逐步收緊。但這會(huì)犧牲解的精確性。分階段求解如果問題可以分解先求解一個(gè)簡(jiǎn)化版如忽略某些非線性約束用其解作為完整問題的初始點(diǎn)。4.2 提供解析梯度與Hessian矩陣默認(rèn)情況下fmincon使用有限差分法來數(shù)值估算目標(biāo)函數(shù)和約束的梯度。這雖然方便但計(jì)算慢且不精確尤其在高維問題中誤差會(huì)放大。顯著提升求解速度和精度的秘訣提供用戶自定義的解析梯度函數(shù)。為目標(biāo)函數(shù)提供梯度創(chuàng)建一個(gè)返回目標(biāo)函數(shù)值f和梯度grad的函數(shù)。function [f, grad] myObjectiveWithGradient(x) f x(1)^2 exp(x(2)); grad [2*x(1); exp(x(2))]; % 梯度向量必須與x同維 end在調(diào)用fmincon時(shí)通過選項(xiàng)啟用并指定梯度函數(shù)options optimoptions(‘fmincon’, ‘SpecifyObjectiveGradient’, true); x fmincon(myObjectiveWithGradient, x0, …, options);為約束提供梯度類似地可以為非線性約束函數(shù)nonlcon提供梯度。這需要函數(shù)返回四個(gè)輸出[c, ceq, gradc, gradceq]其中g(shù)radc和gradceq是約束關(guān)于x的雅可比矩陣轉(zhuǎn)置。設(shè)置options.SpecifyConstraintGradient true。提供Hessian矩陣對(duì)于牛頓類算法提供精確的Hessian矩陣能極大提升收斂速度?？梢酝ㄟ^options.HessianFcn來指定。但對(duì)于擬牛頓法內(nèi)置的BFGS更新已經(jīng)能很好地近似Hessian通常不需要手動(dòng)提供。經(jīng)驗(yàn)之談對(duì)于超過10個(gè)變量的問題強(qiáng)烈建議提供解析梯度。推導(dǎo)梯度雖然需要一些數(shù)學(xué)工作但帶來的性能提升是數(shù)量級(jí)的并且能大大提高求解的魯棒性。使用符號(hào)計(jì)算工具箱Symbolic Math Toolbox可以輔助推導(dǎo)復(fù)雜函數(shù)的梯度。4.3 全局優(yōu)化策略應(yīng)對(duì)非凸難題當(dāng)問題高度非凸時(shí)fmincon只能找到局部最優(yōu)。為了尋找更好的解甚至全局最優(yōu)需要結(jié)合全局優(yōu)化技術(shù)多初始點(diǎn)法這是最簡(jiǎn)單有效的方法。利用循環(huán)或MultiStart對(duì)象從隨機(jī)生成的多個(gè)初始點(diǎn)分別調(diào)用fmincon然后選擇所有結(jié)果中目標(biāo)函數(shù)值最好的那個(gè)。ms MultiStart; problem createOptimProblem(‘fmincon’, ‘objective’, fun, ‘x0’, x0, …); [x_best, fval_best] run(ms, problem, 50); % 從50個(gè)隨機(jī)起點(diǎn)運(yùn)行全局優(yōu)化求解器MATLAB的Global Optimization Toolbox提供了專門的全局優(yōu)化器如ga(遺傳算法)模仿自然選擇適用于變量離散或連續(xù)、問題非光滑的情況。particleswarm(粒子群算法)另一種基于種群的隨機(jī)優(yōu)化方法。simulannealbnd(模擬退火算法)適合變量較少的問題。這些算法通常計(jì)算代價(jià)很高且不能保證找到全局最優(yōu)但能找到比單次局部搜索更好的解。一個(gè)常見的混合策略是先用全局優(yōu)化器如ga進(jìn)行粗略搜索將其結(jié)果作為fmincon的初始點(diǎn)進(jìn)行精細(xì)的局部?jī)?yōu)化。這結(jié)合了全局探索和局部收斂的優(yōu)點(diǎn)。5. 完整案例實(shí)操產(chǎn)品利潤(rùn)最大化模型讓我們通過一個(gè)完整的例子串聯(lián)起所有知識(shí)點(diǎn)。假設(shè)一家公司生產(chǎn)兩種產(chǎn)品其利潤(rùn)函數(shù)單位萬元與產(chǎn)量x1,x2單位千件的關(guān)系為 利潤(rùn) P(x1, x2) 8x1 10x2 - 0.5*(x1^2 x2^2) - 0.2x1x2 生產(chǎn)受到以下限制原材料約束非線性sqrt(x1) 1.5*sqrt(x2) 10機(jī)器工時(shí)約束線性2*x1 3*x2 24市場(chǎng)需求約束線性x1 8, x2 6產(chǎn)量非負(fù)x1 0, x2 0我們的目標(biāo)是最大化利潤(rùn)即最小化-P(x1, x2)。步驟1問題建模與MATLAB代碼實(shí)現(xiàn)% 1. 定義目標(biāo)函數(shù)求最小化所以取負(fù)號(hào) fun (x) -(8*x(1) 10*x(2) - 0.5*(x(1)^2 x(2)^2) - 0.2*x(1)*x(2)); % 2. 定義線性約束 A*x b, Aeq*x beq A [2, 3]; % 2*x1 3*x2 24 b 24; % 無線性等式約束用空數(shù)組表示 Aeq []; beq []; % 3. 定義變量上下界 (lb x ub) lb [0; 0]; ub [8; 6]; % x1 8, x2 6 % 4. 定義非線性約束 sqrt(x1) 1.5*sqrt(x2) 10 function [c, ceq] nonlcon(x) c sqrt(x(1)) 1.5*sqrt(x(2)) - 10; % c 0 ceq []; % 無非線性等式約束 end % 5. 設(shè)置初始猜測(cè)例如中點(diǎn) x0 [4; 3]; % 6. 設(shè)置優(yōu)化選項(xiàng)使用內(nèi)點(diǎn)法顯示迭代過程提高約束容忍度 options optimoptions(‘fmincon’, … ‘Algorithm’, ‘interior-point’, … ‘Display’, ‘iter’, … % 顯示每次迭代信息 ‘ConstraintTolerance’, 1e-8, … % 約束容忍度 ‘OptimalityTolerance’, 1e-8); % 最優(yōu)性容忍度 % 7. 調(diào)用 fmincon 求解 [x_opt, fval_opt, exitflag, output] fmincon(fun, x0, A, b, Aeq, beq, lb, ub, nonlcon, options); % 8. 輸出結(jié)果 fprintf(‘最優(yōu)產(chǎn)量x1 %.4f (千件), x2 %.4f (千件)\n’, x_opt(1), x_opt(2)); fprintf(‘最大利潤(rùn)%.4f (萬元)\n’, -fval_opt); % 注意取負(fù)號(hào)轉(zhuǎn)回利潤(rùn) fprintf(‘退出標(biāo)志%d\n’, exitflag); fprintf(‘迭代次數(shù)%d\n’, output.iterations); fprintf(‘一階最優(yōu)性度量%.2e\n’, output.firstorderopt); % 9. 驗(yàn)證約束 fprintf(‘\n約束驗(yàn)證\n’); fprintf(‘原材料約束sqrt(x1)1.5*sqrt(x2) %.4f 10\n’, sqrt(x_opt(1)) 1.5*sqrt(x_opt(2))); fprintf(‘機(jī)器工時(shí)約束2*x13*x2 %.4f 24\n’, 2*x_opt(1)3*x_opt(2)); fprintf(‘市場(chǎng)需求約束x1%.4f8, x2%.4f6\n’, x_opt(1), x_opt(2));步驟2結(jié)果分析與解讀運(yùn)行上述代碼你會(huì)看到類似以下的迭代輸出和結(jié)果Iter F-count f(x) Feasibility Steplength Step First-order optimality 0 3 -3.220000e01 1.000e00 1 6 -3.496263e01 0.000e00 1.000e00 1.604e00 1.053e00 2 9 -3.496263e01 0.000e00 1.000e00 1.604e00 1.053e00 ... 最優(yōu)產(chǎn)量x1 4.0000 (千件), x2 5.3333 (千件) 最大利潤(rùn)34.9626 (萬元) 退出標(biāo)志1 迭代次數(shù)8 一階最優(yōu)性度量1.05e-06 約束驗(yàn)證 原材料約束sqrt(x1)1.5*sqrt(x2) 10.0000 10 緊約束 機(jī)器工時(shí)約束2*x13*x2 24.0000 24 緊約束 市場(chǎng)需求約束x14.00008, x25.33336分析退出標(biāo)志為1說明算法成功收斂到一個(gè)局部最優(yōu)解對(duì)于此凸問題也是全局最優(yōu)。兩個(gè)約束原材料和機(jī)器工時(shí)在最優(yōu)解處都是“緊”的等號(hào)成立這意味著這些資源被完全利用是限制利潤(rùn)增長(zhǎng)的關(guān)鍵瓶頸。市場(chǎng)需求約束并未達(dá)到上限說明不是當(dāng)前生產(chǎn)計(jì)劃的限制因素。一階最優(yōu)性度量非常小1.05e-06遠(yuǎn)小于默認(rèn)容差1e-6說明解的質(zhì)量很高。從迭代過程看算法很快找到了可行域Feasibility從1變?yōu)?并在幾步內(nèi)收斂。步驟3敏感性分析與“What-If”建模的價(jià)值不止于得到一個(gè)數(shù)字。我們可以利用這個(gè)模型進(jìn)行簡(jiǎn)單的敏感性分析如果原材料供應(yīng)增加10%會(huì)怎樣將非線性約束的右端項(xiàng)從10改為11重新求解。你會(huì)發(fā)現(xiàn)利潤(rùn)增加了并且可能某個(gè)之前“緊”的約束變得“松”了這能指導(dǎo)采購(gòu)決策。如果產(chǎn)品2的市場(chǎng)需求上限提高到7呢修改ub(2) 7重新求解。觀察最優(yōu)產(chǎn)量x2是否增加以及利潤(rùn)的提升幅度這能評(píng)估市場(chǎng)擴(kuò)張的潛在收益。踩坑記錄在這個(gè)例子中非線性約束涉及sqrt(x)。必須確保初始點(diǎn)x0和迭代過程中的x不會(huì)為負(fù)否則sqrt會(huì)返回復(fù)數(shù)導(dǎo)致求解失敗。這就是為什么我們?cè)O(shè)置了lb [0; 0]。在實(shí)際問題中遇到對(duì)數(shù)函數(shù)log(x)、分?jǐn)?shù)冪等同樣要特別注意定義域通過設(shè)置合理的下界來保證數(shù)值穩(wěn)定性。6. 常見問題排查與調(diào)試技巧即使按照指南操作在實(shí)際編碼和求解中依然會(huì)遇到各種問題。這里匯總了一些典型錯(cuò)誤及其解決方法。6.1 求解失敗或結(jié)果不理想問題現(xiàn)象可能原因排查與解決步驟exitflag為負(fù)數(shù)1. 目標(biāo)函數(shù)或約束函數(shù)返回NaN/Inf。2. 初始點(diǎn)x0不可行對(duì)某些算法。3. 問題可能無界。1.添加調(diào)試輸出在自定義函數(shù)開頭添加disp(x)或在函數(shù)內(nèi)設(shè)置斷點(diǎn)檢查導(dǎo)致非數(shù)值的輸入。2.檢查定義域確保log,sqrt, 除法等運(yùn)算在定義域內(nèi)。3.嘗試一個(gè)更可行的初始點(diǎn)或使用‘interior-point’算法它對(duì)初始可行性要求較低。4. 檢查模型邏輯目標(biāo)函數(shù)是否可能無限減小。exitflag為 0迭代次數(shù)或函數(shù)計(jì)算次數(shù)達(dá)到上限。1. 增加options.MaxIterations和options.MaxFunctionEvaluations。2. 檢查是否收斂緩慢。提供解析梯度通常能極大加速收斂。3. 嘗試不同的算法如從‘interior-point’切換到‘sqp’。解不滿足約束約束容忍度 (ConstraintTolerance) 設(shè)置得過大或者算法在數(shù)值誤差下提前終止。1. 檢查output.constrviolation查看最大約束違反值。2. 減小options.ConstraintTolerance例如1e-8。3. 手動(dòng)驗(yàn)證解是否滿足約束如案例中所做。每次運(yùn)行結(jié)果差異大問題是非凸的算法收斂到不同的局部最優(yōu)解。1. 使用MultiStart從多個(gè)隨機(jī)初始點(diǎn)求解。2. 考慮使用全局優(yōu)化算法如ga進(jìn)行初步搜索。求解速度極慢1. 目標(biāo)函數(shù)/約束函數(shù)本身計(jì)算復(fù)雜。2. 使用有限差分計(jì)算梯度高維問題尤甚。3. 問題規(guī)模太大。1.優(yōu)化函數(shù)代碼向量化操作避免循環(huán)。2.提供解析梯度這是提升速度最有效的方法。3. 對(duì)于大規(guī)模問題確保使用‘interior-point’算法并利用稀疏矩陣。6.2 數(shù)值穩(wěn)定性與技巧縮放變量如果決策變量的數(shù)量級(jí)相差巨大如x1約1e-6,x2約1e3會(huì)導(dǎo)致Hessian矩陣條件數(shù)很差嚴(yán)重影響算法數(shù)值穩(wěn)定性。最佳實(shí)踐是對(duì)變量進(jìn)行縮放使其數(shù)量級(jí)大致在1附近。例如定義新變量y1 1e6 * x1,y2 1e-3 * x2在模型中用y代替x求解后再轉(zhuǎn)換回來。避免數(shù)值微分如前所述盡量提供解析導(dǎo)數(shù)。如果實(shí)在無法推導(dǎo)可以考慮使用自動(dòng)微分AD工具但對(duì)于MATLAB用戶提供解析梯度是最直接的。檢查梯度在提供自定義梯度后務(wù)必使用options.CheckGradients true進(jìn)行驗(yàn)證。MATLAB會(huì)將你的梯度函數(shù)計(jì)算結(jié)果與有限差分結(jié)果進(jìn)行比較并報(bào)告差異。這是排除梯度計(jì)算錯(cuò)誤的關(guān)鍵一步。理解“容差”O(jiān)ptimalityTolerance和StepTolerance決定了算法何時(shí)停止。通常1e-6是默認(rèn)且合理的值。對(duì)于工程應(yīng)用1e-4可能已足夠精確。過分追求1e-12這樣的高精度只會(huì)無謂增加計(jì)算時(shí)間。非線性規(guī)劃是連接數(shù)學(xué)模型與現(xiàn)實(shí)復(fù)雜決策的橋梁。它沒有銀彈其魅力在于需要你根據(jù)具體問題的“脾氣”凸性、光滑性、規(guī)模來選擇合適的“工具”和“策略”。從理解問題分類開始到熟練運(yùn)用fmincon的各項(xiàng)功能再到掌握提供梯度、處理非凸、調(diào)試錯(cuò)誤等進(jìn)階技巧每一步都伴隨著從理論到實(shí)踐的深化。記住一個(gè)成功的求解往往始于一個(gè)合理的模型表述依賴于一個(gè)明智的算法和選項(xiàng)配置并最終得益于對(duì)結(jié)果的嚴(yán)謹(jǐn)驗(yàn)證和敏感性分析。多動(dòng)手、多試錯(cuò)、多思考“為什么”你就能將這門技術(shù)真正化為解決實(shí)際難題的利器。 MATLAB非線性規(guī)劃實(shí)戰(zhàn)：從fmincon算法選擇到全局優(yōu)化策略
1. 項(xiàng)目概述從線性到非線性的思維躍遷在數(shù)學(xué)建模的實(shí)戰(zhàn)中我們遇到的絕大多數(shù)問題其目標(biāo)函數(shù)或約束條件都不是簡(jiǎn)單的線性關(guān)系。比如你想優(yōu)化一個(gè)工廠的生產(chǎn)計(jì)劃成本可能隨著產(chǎn)量呈指數(shù)增長(zhǎng)目標(biāo)函數(shù)非線性或者你設(shè)計(jì)一個(gè)機(jī)械結(jié)構(gòu)其應(yīng)力必須小于材料的非線性屈服強(qiáng)度約束條件非線性。這時(shí)線性規(guī)劃那套漂亮的單純形法就完全失效了。非線性規(guī)劃正是為了解決這類“彎彎繞繞”的優(yōu)化問題而生的核心數(shù)學(xué)工具。它不像線性規(guī)劃那樣有“標(biāo)準(zhǔn)答案”式的通用解法更像是一個(gè)工具箱里面裝著各種針對(duì)不同問題特性的“專用扳手”。我接觸過很多剛開始做建模的同學(xué)一看到“非線性”三個(gè)字就頭疼覺得深不可測(cè)。其實(shí)不然它的核心思想非常直觀在復(fù)雜的地形目標(biāo)函數(shù)曲面上找到那個(gè)最低點(diǎn)最小值或最高點(diǎn)最大值同時(shí)不能跑到禁區(qū)約束條件里去。這次我們就來徹底拆解這個(gè)工具箱不僅告訴你每個(gè)工具算法怎么用更重點(diǎn)講清楚什么時(shí)候該用哪個(gè)以及用的時(shí)候最容易在哪兒翻車。我們會(huì)以最常用的MATLAB環(huán)境為例手把手帶你從理論走到代碼實(shí)現(xiàn)讓你下次遇到非線性問題時(shí)能胸有成竹地選出最合適的那把“扳手”。2. 非線性規(guī)劃的核心思想與問題分類在動(dòng)手寫代碼之前我們必須先搞清楚面對(duì)的是什么“型號(hào)”的問題。非線性規(guī)劃問題通常可以寫成如下標(biāo)準(zhǔn)形式最小化問題Minimize: f(x) Subject to: g_i(x) ≤ 0, i 1, ..., m (不等式約束) h_j(x) 0, j 1, ..., p (等式約束) x ∈ R^n (決策變量)這里f(x), g_i(x), h_j(x) 中至少有一個(gè)是非線性函數(shù)。根據(jù)這些函數(shù)的特性我們可以把問題分門別類這直接決定了我們?cè)撨x用哪種算法。2.1 凸與非凸決定問題難度的分水嶺這是非線性規(guī)劃中最關(guān)鍵的分類沒有之一。凸規(guī)劃如果目標(biāo)函數(shù) f(x) 是凸函數(shù)并且不等式約束函數(shù) g_i(x) 是凸函數(shù)等式約束 h_j(x) 是線性函數(shù)那么這個(gè)問題就是凸規(guī)劃。凸規(guī)劃的任何局部最優(yōu)解必定是全局最優(yōu)解。這是它最大的優(yōu)點(diǎn)意味著算法只要找到一個(gè)“坑底”那就是整個(gè)區(qū)域的最低點(diǎn)。求解凸規(guī)劃相對(duì)“友好”。非凸規(guī)劃不滿足上述凸性條件的規(guī)劃問題。它的“地形圖”可能像連綿的群山有無數(shù)個(gè)山谷局部最優(yōu)點(diǎn)算法很容易陷在某個(gè)小山谷里而找不到最深的那一個(gè)全局最優(yōu)點(diǎn)。求解非凸規(guī)劃是NP-Hard問題通常只能尋找“較好的”局部最優(yōu)解或者采用一些隨機(jī)策略如模擬退火、遺傳算法來嘗試尋找全局最優(yōu)。實(shí)操心得在實(shí)際建模中我們首先應(yīng)該嘗試判斷問題是否具有凸性。一個(gè)簡(jiǎn)單的技巧如果目標(biāo)函數(shù)是二次型且Hessian矩陣半正定或者約束是線性的那么它很可能是凸的。對(duì)于復(fù)雜函數(shù)判斷凸性需要利用二階條件Hessian矩陣處處半正定這在實(shí)踐中往往很困難。因此一個(gè)務(wù)實(shí)的做法是默認(rèn)問題是非凸的然后選擇能處理非凸問題的穩(wěn)健算法同時(shí)嘗試從多個(gè)不同的初始點(diǎn)出發(fā)求解以降低陷入糟糕局部最優(yōu)的風(fēng)險(xiǎn)。2.2 無約束與有約束解決問題的基本框架無約束非線性優(yōu)化問題中沒有任何 g_i(x) 和 h_j(x) 的限制。這類問題的經(jīng)典算法構(gòu)成了非線性優(yōu)化的基石例如梯度下降法沿著目標(biāo)函數(shù)負(fù)梯度方向迭代簡(jiǎn)單但收斂慢。牛頓法利用目標(biāo)函數(shù)的二階導(dǎo)數(shù)Hessian矩陣信息收斂速度快但需要計(jì)算Hessian矩陣及其逆計(jì)算量大。擬牛頓法如BFGS, DFP通過構(gòu)造一個(gè)近似矩陣來模擬Hessian矩陣的逆既保持了較快的收斂速度又避免了直接計(jì)算Hessian矩陣是實(shí)踐中無約束優(yōu)化的首選。有約束非線性優(yōu)化這是我們討論的重點(diǎn)也是fmincon等求解器主要應(yīng)對(duì)的場(chǎng)景。核心思路是將有約束問題轉(zhuǎn)化為一系列無約束或更簡(jiǎn)單的約束問題來求解。2.3 二次規(guī)劃非線性中的“線性”特例二次規(guī)劃是指目標(biāo)函數(shù)是二次函數(shù)約束條件是線性函數(shù)的一類特殊非線性規(guī)劃。它的標(biāo)準(zhǔn)形式為 Minimize: (1/2) * x^T * H * x c^T * x Subject to: A * x ≤ b, Aeq * x beq, lb ≤ x ≤ ub雖然目標(biāo)函數(shù)是非線性的二次但由于其結(jié)構(gòu)特殊存在非常高效和可靠的專用算法如有效集法、內(nèi)點(diǎn)法。在MATLAB中可以使用quadprog函數(shù)專門求解QP問題。很多復(fù)雜的非線性問題在局部可以用二次函數(shù)來近似因此QP求解器也是許多高級(jí)非線性算法如序列二次規(guī)劃SQP的核心子步驟。3. MATLAB實(shí)戰(zhàn)核心fmincon求解器深度解析MATLAB的fmincon是求解中小規(guī)模有約束非線性規(guī)劃問題的“瑞士軍刀”。它的強(qiáng)大之處在于內(nèi)部集成了多種算法可以自動(dòng)或手動(dòng)選擇以適應(yīng)不同問題。但要用好它必須理解其每一個(gè)參數(shù)背后的意義。3.1fmincon的基本調(diào)用與參數(shù)精講一個(gè)最基礎(chǔ)的調(diào)用格式如下[x, fval, exitflag, output] fmincon(fun, x0, A, b, Aeq, beq, lb, ub, nonlcon, options)我們來逐一拆解這些輸入輸出參數(shù)fun目標(biāo)函數(shù)句柄。例如(x) x(1)^2 x(2)^2。這里最容易出錯(cuò)的地方是函數(shù)定義必須能接受向量輸入x并返回標(biāo)量值。如果目標(biāo)函數(shù)計(jì)算量很大可以考慮在函數(shù)內(nèi)部進(jìn)行向量化操作或使用全局變量/嵌套函數(shù)傳遞額外參數(shù)。x0初始猜測(cè)值。這是影響求解結(jié)果最關(guān)鍵的因素之一尤其對(duì)于非凸問題。糟糕的初始點(diǎn)可能導(dǎo)致算法收斂到很差的局部最優(yōu)甚至失敗。一個(gè)好的策略是根據(jù)物理意義或經(jīng)驗(yàn)給出初始值或者進(jìn)行簡(jiǎn)單的網(wǎng)格搜索從多個(gè)初始點(diǎn)中選取最好的結(jié)果。A, b, Aeq, beq, lb, ub線性約束和邊界約束。這是定義約束最高效的方式應(yīng)優(yōu)先使用。nonlcon非線性約束函數(shù)句柄。該函數(shù)需要返回兩個(gè)輸出[c, ceq]其中c(x) 0表示非線性不等式約束ceq(x) 0表示非線性等式約束。即使只有一種約束也必須同時(shí)返回兩個(gè)輸出將不存在的那個(gè)設(shè)為空數(shù)組[]。3.2 算法選擇options的設(shè)置藝術(shù)通過optimoptions(‘fmincon’)來設(shè)置選項(xiàng)。算法選擇 (Algorithm) 是核心‘interior-point’(內(nèi)點(diǎn)法)默認(rèn)且最通用的算法。特別適合大規(guī)模問題能高效處理邊界約束和稀疏性。它通過在可行域內(nèi)部構(gòu)造一條路徑逼近最優(yōu)解。對(duì)于大多數(shù)問題首選這個(gè)算法?！畇qp’(序列二次規(guī)劃)另一種強(qiáng)大的通用算法。它通過在每一步迭代中求解一個(gè)二次規(guī)劃子問題來尋找搜索方向。對(duì)于中小規(guī)模問題尤其是約束較多的問題表現(xiàn)可能比內(nèi)點(diǎn)法更好?！產(chǎn)ctive-set’(有效集法)一種較老的算法適用于中小規(guī)模問題。它能精確識(shí)別在最優(yōu)解處起作用的約束active constraints。對(duì)于需要知道哪些約束是“緊”的問題這個(gè)算法有優(yōu)勢(shì)?！畉rust-region-reflective’(信賴域反射法)這個(gè)算法要求目標(biāo)函數(shù)是非線性的標(biāo)量函數(shù)且只能處理邊界約束或線性等式約束不能直接處理非線性約束或線性不等式約束但可以通過轉(zhuǎn)換。它的優(yōu)勢(shì)在于能利用目標(biāo)函數(shù)的梯度信息對(duì)于特定類型的問題非常高效。避坑指南如果你的問題包含非線性約束那么算法只能從‘interior-point’和‘sqp’中選擇。初次求解一個(gè)未知問題時(shí)建議先使用默認(rèn)的‘interior-point’算法。如果收斂速度慢或不穩(wěn)定再嘗試‘sqp’。務(wù)必在選項(xiàng)中打開梯度檢查options optimoptions(‘fmincon’, ‘CheckGradients’, true);這能幫你發(fā)現(xiàn)自定義梯度函數(shù)中的錯(cuò)誤避免因梯度不準(zhǔn)導(dǎo)致算法失敗。3.3 輸出結(jié)果解讀與診斷求解完成后不能只看最優(yōu)解x和最優(yōu)值fvalexitflag和output包含了至關(guān)重要的診斷信息。exitflag(退出標(biāo)志) 0算法收斂到局部最優(yōu)解。這是成功標(biāo)志。 0迭代次數(shù)或函數(shù)計(jì)算次數(shù)超過了MaxIterations或MaxFunctionEvaluations選項(xiàng)設(shè)置的最大值。此時(shí)得到的解可能不是最優(yōu)的需要增加迭代上限或檢查問題 formulation。 0求解失敗。常見原因包括目標(biāo)函數(shù)或約束函數(shù)在迭代點(diǎn)處返回了NaN或Inf問題可能無界初始點(diǎn)不可行對(duì)于嚴(yán)格要求可行性的算法。需要根據(jù)output.message中的信息進(jìn)行排查。output結(jié)構(gòu)體包含迭代次數(shù) (iterations)、函數(shù)計(jì)算次數(shù) (funcCount)、一階最優(yōu)性條件 (firstorderopt)、算法類型等。firstorderopt衡量了當(dāng)前解滿足一階最優(yōu)性條件KKT條件的程度這個(gè)值越小說明解越“優(yōu)”。通常小于1e-6可以認(rèn)為是很好的收斂。4. 進(jìn)階技術(shù)與實(shí)戰(zhàn)策略掌握了fmincon的基本用法我們來看看如何解決更復(fù)雜的情況以及如何提升求解的效率和穩(wěn)定性。4.1 處理復(fù)雜非線性約束與可行性當(dāng)非線性約束非常復(fù)雜時(shí)算法可能很難找到一個(gè)可行的初始點(diǎn)或者在迭代中保持可行性。這時(shí)可以嘗試使用罰函數(shù)法這是將約束問題轉(zhuǎn)化為無約束問題的經(jīng)典思路?；舅枷胧菍⒓s束違反的程度作為一個(gè)“懲罰項(xiàng)”加到目標(biāo)函數(shù)中。例如對(duì)于約束g(x) 0可以構(gòu)造罰函數(shù)P(x) f(x) μ * max(0, g(x))^2其中μ是一個(gè)很大的正數(shù)罰因子。然后使用無約束優(yōu)化方法如fminunc求解P(x)。隨著μ增大解會(huì)越來越逼近原約束問題的解。缺點(diǎn)是罰因子需要精心選擇太大可能導(dǎo)致數(shù)值問題太小則約束得不到滿足。fmincon的可行性模式對(duì)于某些算法如sqp可以設(shè)置options.ConstraintTolerance來放寬對(duì)約束的嚴(yán)格滿足要求讓算法先找到一個(gè)“差不多”可行的點(diǎn)再逐步收緊。但這會(huì)犧牲解的精確性。分階段求解如果問題可以分解先求解一個(gè)簡(jiǎn)化版如忽略某些非線性約束用其解作為完整問題的初始點(diǎn)。4.2 提供解析梯度與Hessian矩陣默認(rèn)情況下fmincon使用有限差分法來數(shù)值估算目標(biāo)函數(shù)和約束的梯度。這雖然方便但計(jì)算慢且不精確尤其在高維問題中誤差會(huì)放大。顯著提升求解速度和精度的秘訣提供用戶自定義的解析梯度函數(shù)。為目標(biāo)函數(shù)提供梯度創(chuàng)建一個(gè)返回目標(biāo)函數(shù)值f和梯度grad的函數(shù)。function [f, grad] myObjectiveWithGradient(x) f x(1)^2 exp(x(2)); grad [2*x(1); exp(x(2))]; % 梯度向量必須與x同維 end在調(diào)用fmincon時(shí)通過選項(xiàng)啟用并指定梯度函數(shù)options optimoptions(‘fmincon’, ‘SpecifyObjectiveGradient’, true); x fmincon(myObjectiveWithGradient, x0, …, options);為約束提供梯度類似地可以為非線性約束函數(shù)nonlcon提供梯度。這需要函數(shù)返回四個(gè)輸出[c, ceq, gradc, gradceq]其中g(shù)radc和gradceq是約束關(guān)于x的雅可比矩陣轉(zhuǎn)置。設(shè)置options.SpecifyConstraintGradient true。提供Hessian矩陣對(duì)于牛頓類算法提供精確的Hessian矩陣能極大提升收斂速度?？梢酝ㄟ^options.HessianFcn來指定。但對(duì)于擬牛頓法內(nèi)置的BFGS更新已經(jīng)能很好地近似Hessian通常不需要手動(dòng)提供。經(jīng)驗(yàn)之談對(duì)于超過10個(gè)變量的問題強(qiáng)烈建議提供解析梯度。推導(dǎo)梯度雖然需要一些數(shù)學(xué)工作但帶來的性能提升是數(shù)量級(jí)的并且能大大提高求解的魯棒性。使用符號(hào)計(jì)算工具箱Symbolic Math Toolbox可以輔助推導(dǎo)復(fù)雜函數(shù)的梯度。4.3 全局優(yōu)化策略應(yīng)對(duì)非凸難題當(dāng)問題高度非凸時(shí)fmincon只能找到局部最優(yōu)。為了尋找更好的解甚至全局最優(yōu)需要結(jié)合全局優(yōu)化技術(shù)多初始點(diǎn)法這是最簡(jiǎn)單有效的方法。利用循環(huán)或MultiStart對(duì)象從隨機(jī)生成的多個(gè)初始點(diǎn)分別調(diào)用fmincon然后選擇所有結(jié)果中目標(biāo)函數(shù)值最好的那個(gè)。ms MultiStart; problem createOptimProblem(‘fmincon’, ‘objective’, fun, ‘x0’, x0, …); [x_best, fval_best] run(ms, problem, 50); % 從50個(gè)隨機(jī)起點(diǎn)運(yùn)行全局優(yōu)化求解器MATLAB的Global Optimization Toolbox提供了專門的全局優(yōu)化器如ga(遺傳算法)模仿自然選擇適用于變量離散或連續(xù)、問題非光滑的情況。particleswarm(粒子群算法)另一種基于種群的隨機(jī)優(yōu)化方法。simulannealbnd(模擬退火算法)適合變量較少的問題。這些算法通常計(jì)算代價(jià)很高且不能保證找到全局最優(yōu)但能找到比單次局部搜索更好的解。一個(gè)常見的混合策略是先用全局優(yōu)化器如ga進(jìn)行粗略搜索將其結(jié)果作為fmincon的初始點(diǎn)進(jìn)行精細(xì)的局部?jī)?yōu)化。這結(jié)合了全局探索和局部收斂的優(yōu)點(diǎn)。5. 完整案例實(shí)操產(chǎn)品利潤(rùn)最大化模型讓我們通過一個(gè)完整的例子串聯(lián)起所有知識(shí)點(diǎn)。假設(shè)一家公司生產(chǎn)兩種產(chǎn)品其利潤(rùn)函數(shù)單位萬元與產(chǎn)量x1,x2單位千件的關(guān)系為 利潤(rùn) P(x1, x2) 8x1 10x2 - 0.5*(x1^2 x2^2) - 0.2x1x2 生產(chǎn)受到以下限制原材料約束非線性sqrt(x1) 1.5*sqrt(x2) 10機(jī)器工時(shí)約束線性2*x1 3*x2 24市場(chǎng)需求約束線性x1 8, x2 6產(chǎn)量非負(fù)x1 0, x2 0我們的目標(biāo)是最大化利潤(rùn)即最小化-P(x1, x2)。步驟1問題建模與MATLAB代碼實(shí)現(xiàn)% 1. 定義目標(biāo)函數(shù)求最小化所以取負(fù)號(hào) fun (x) -(8*x(1) 10*x(2) - 0.5*(x(1)^2 x(2)^2) - 0.2*x(1)*x(2)); % 2. 定義線性約束 A*x b, Aeq*x beq A [2, 3]; % 2*x1 3*x2 24 b 24; % 無線性等式約束用空數(shù)組表示 Aeq []; beq []; % 3. 定義變量上下界 (lb x ub) lb [0; 0]; ub [8; 6]; % x1 8, x2 6 % 4. 定義非線性約束 sqrt(x1) 1.5*sqrt(x2) 10 function [c, ceq] nonlcon(x) c sqrt(x(1)) 1.5*sqrt(x(2)) - 10; % c 0 ceq []; % 無非線性等式約束 end % 5. 設(shè)置初始猜測(cè)例如中點(diǎn) x0 [4; 3]; % 6. 設(shè)置優(yōu)化選項(xiàng)使用內(nèi)點(diǎn)法顯示迭代過程提高約束容忍度 options optimoptions(‘fmincon’, … ‘Algorithm’, ‘interior-point’, … ‘Display’, ‘iter’, … % 顯示每次迭代信息 ‘ConstraintTolerance’, 1e-8, … % 約束容忍度 ‘OptimalityTolerance’, 1e-8); % 最優(yōu)性容忍度 % 7. 調(diào)用 fmincon 求解 [x_opt, fval_opt, exitflag, output] fmincon(fun, x0, A, b, Aeq, beq, lb, ub, nonlcon, options); % 8. 輸出結(jié)果 fprintf(‘最優(yōu)產(chǎn)量x1 %.4f (千件), x2 %.4f (千件)\n’, x_opt(1), x_opt(2)); fprintf(‘最大利潤(rùn)%.4f (萬元)\n’, -fval_opt); % 注意取負(fù)號(hào)轉(zhuǎn)回利潤(rùn) fprintf(‘退出標(biāo)志%d\n’, exitflag); fprintf(‘迭代次數(shù)%d\n’, output.iterations); fprintf(‘一階最優(yōu)性度量%.2e\n’, output.firstorderopt); % 9. 驗(yàn)證約束 fprintf(‘\n約束驗(yàn)證\n’); fprintf(‘原材料約束sqrt(x1)1.5*sqrt(x2) %.4f 10\n’, sqrt(x_opt(1)) 1.5*sqrt(x_opt(2))); fprintf(‘機(jī)器工時(shí)約束2*x13*x2 %.4f 24\n’, 2*x_opt(1)3*x_opt(2)); fprintf(‘市場(chǎng)需求約束x1%.4f8, x2%.4f6\n’, x_opt(1), x_opt(2));步驟2結(jié)果分析與解讀運(yùn)行上述代碼你會(huì)看到類似以下的迭代輸出和結(jié)果Iter F-count f(x) Feasibility Steplength Step First-order optimality 0 3 -3.220000e01 1.000e00 1 6 -3.496263e01 0.000e00 1.000e00 1.604e00 1.053e00 2 9 -3.496263e01 0.000e00 1.000e00 1.604e00 1.053e00 ... 最優(yōu)產(chǎn)量x1 4.0000 (千件), x2 5.3333 (千件) 最大利潤(rùn)34.9626 (萬元) 退出標(biāo)志1 迭代次數(shù)8 一階最優(yōu)性度量1.05e-06 約束驗(yàn)證 原材料約束sqrt(x1)1.5*sqrt(x2) 10.0000 10 緊約束 機(jī)器工時(shí)約束2*x13*x2 24.0000 24 緊約束 市場(chǎng)需求約束x14.00008, x25.33336分析退出標(biāo)志為1說明算法成功收斂到一個(gè)局部最優(yōu)解對(duì)于此凸問題也是全局最優(yōu)。兩個(gè)約束原材料和機(jī)器工時(shí)在最優(yōu)解處都是“緊”的等號(hào)成立這意味著這些資源被完全利用是限制利潤(rùn)增長(zhǎng)的關(guān)鍵瓶頸。市場(chǎng)需求約束并未達(dá)到上限說明不是當(dāng)前生產(chǎn)計(jì)劃的限制因素。一階最優(yōu)性度量非常小1.05e-06遠(yuǎn)小于默認(rèn)容差1e-6說明解的質(zhì)量很高。從迭代過程看算法很快找到了可行域Feasibility從1變?yōu)?并在幾步內(nèi)收斂。步驟3敏感性分析與“What-If”建模的價(jià)值不止于得到一個(gè)數(shù)字。我們可以利用這個(gè)模型進(jìn)行簡(jiǎn)單的敏感性分析如果原材料供應(yīng)增加10%會(huì)怎樣將非線性約束的右端項(xiàng)從10改為11重新求解。你會(huì)發(fā)現(xiàn)利潤(rùn)增加了并且可能某個(gè)之前“緊”的約束變得“松”了這能指導(dǎo)采購(gòu)決策。如果產(chǎn)品2的市場(chǎng)需求上限提高到7呢修改ub(2) 7重新求解。觀察最優(yōu)產(chǎn)量x2是否增加以及利潤(rùn)的提升幅度這能評(píng)估市場(chǎng)擴(kuò)張的潛在收益。踩坑記錄在這個(gè)例子中非線性約束涉及sqrt(x)。必須確保初始點(diǎn)x0和迭代過程中的x不會(huì)為負(fù)否則sqrt會(huì)返回復(fù)數(shù)導(dǎo)致求解失敗。這就是為什么我們?cè)O(shè)置了lb [0; 0]。在實(shí)際問題中遇到對(duì)數(shù)函數(shù)log(x)、分?jǐn)?shù)冪等同樣要特別注意定義域通過設(shè)置合理的下界來保證數(shù)值穩(wěn)定性。6. 常見問題排查與調(diào)試技巧即使按照指南操作在實(shí)際編碼和求解中依然會(huì)遇到各種問題。這里匯總了一些典型錯(cuò)誤及其解決方法。6.1 求解失敗或結(jié)果不理想問題現(xiàn)象可能原因排查與解決步驟exitflag為負(fù)數(shù)1. 目標(biāo)函數(shù)或約束函數(shù)返回NaN/Inf。2. 初始點(diǎn)x0不可行對(duì)某些算法。3. 問題可能無界。1.添加調(diào)試輸出在自定義函數(shù)開頭添加disp(x)或在函數(shù)內(nèi)設(shè)置斷點(diǎn)檢查導(dǎo)致非數(shù)值的輸入。2.檢查定義域確保log,sqrt, 除法等運(yùn)算在定義域內(nèi)。3.嘗試一個(gè)更可行的初始點(diǎn)或使用‘interior-point’算法它對(duì)初始可行性要求較低。4. 檢查模型邏輯目標(biāo)函數(shù)是否可能無限減小。exitflag為 0迭代次數(shù)或函數(shù)計(jì)算次數(shù)達(dá)到上限。1. 增加options.MaxIterations和options.MaxFunctionEvaluations。2. 檢查是否收斂緩慢。提供解析梯度通常能極大加速收斂。3. 嘗試不同的算法如從‘interior-point’切換到‘sqp’。解不滿足約束約束容忍度 (ConstraintTolerance) 設(shè)置得過大或者算法在數(shù)值誤差下提前終止。1. 檢查output.constrviolation查看最大約束違反值。2. 減小options.ConstraintTolerance例如1e-8。3. 手動(dòng)驗(yàn)證解是否滿足約束如案例中所做。每次運(yùn)行結(jié)果差異大問題是非凸的算法收斂到不同的局部最優(yōu)解。1. 使用MultiStart從多個(gè)隨機(jī)初始點(diǎn)求解。2. 考慮使用全局優(yōu)化算法如ga進(jìn)行初步搜索。求解速度極慢1. 目標(biāo)函數(shù)/約束函數(shù)本身計(jì)算復(fù)雜。2. 使用有限差分計(jì)算梯度高維問題尤甚。3. 問題規(guī)模太大。1.優(yōu)化函數(shù)代碼向量化操作避免循環(huán)。2.提供解析梯度這是提升速度最有效的方法。3. 對(duì)于大規(guī)模問題確保使用‘interior-point’算法并利用稀疏矩陣。6.2 數(shù)值穩(wěn)定性與技巧縮放變量如果決策變量的數(shù)量級(jí)相差巨大如x1約1e-6,x2約1e3會(huì)導(dǎo)致Hessian矩陣條件數(shù)很差嚴(yán)重影響算法數(shù)值穩(wěn)定性。最佳實(shí)踐是對(duì)變量進(jìn)行縮放使其數(shù)量級(jí)大致在1附近。例如定義新變量y1 1e6 * x1,y2 1e-3 * x2在模型中用y代替x求解后再轉(zhuǎn)換回來。避免數(shù)值微分如前所述盡量提供解析導(dǎo)數(shù)。如果實(shí)在無法推導(dǎo)可以考慮使用自動(dòng)微分AD工具但對(duì)于MATLAB用戶提供解析梯度是最直接的。檢查梯度在提供自定義梯度后務(wù)必使用options.CheckGradients true進(jìn)行驗(yàn)證。MATLAB會(huì)將你的梯度函數(shù)計(jì)算結(jié)果與有限差分結(jié)果進(jìn)行比較并報(bào)告差異。這是排除梯度計(jì)算錯(cuò)誤的關(guān)鍵一步。理解“容差”O(jiān)ptimalityTolerance和StepTolerance決定了算法何時(shí)停止。通常1e-6是默認(rèn)且合理的值。對(duì)于工程應(yīng)用1e-4可能已足夠精確。過分追求1e-12這樣的高精度只會(huì)無謂增加計(jì)算時(shí)間。非線性規(guī)劃是連接數(shù)學(xué)模型與現(xiàn)實(shí)復(fù)雜決策的橋梁。它沒有銀彈其魅力在于需要你根據(jù)具體問題的“脾氣”凸性、光滑性、規(guī)模來選擇合適的“工具”和“策略”。從理解問題分類開始到熟練運(yùn)用fmincon的各項(xiàng)功能再到掌握提供梯度、處理非凸、調(diào)試錯(cuò)誤等進(jìn)階技巧每一步都伴隨著從理論到實(shí)踐的深化。記住一個(gè)成功的求解往往始于一個(gè)合理的模型表述依賴于一個(gè)明智的算法和選項(xiàng)配置并最終得益于對(duì)結(jié)果的嚴(yán)謹(jǐn)驗(yàn)證和敏感性分析。多動(dòng)手、多試錯(cuò)、多思考“為什么”你就能將這門技術(shù)真正化為解決實(shí)際難題的利器。