第五講Matlab優(yōu)化工具箱簡(jiǎn)介 optimizationtoobox 1 線性優(yōu)化2 非線性優(yōu)化3 極小化極大 Minmax 問(wèn)題4 曲線擬合與插值 線性規(guī)劃問(wèn)題是目標(biāo)函數(shù)和約束條件均為線性函數(shù)的問(wèn)題 MATLAB7 0解決的線性規(guī)劃問(wèn)題的標(biāo)準(zhǔn)形式為minsub to 其中f x b beq lb ub為向量 A Aeq為矩陣 其它形式的線性規(guī)劃問(wèn)題都可經(jīng)過(guò)適當(dāng)變換化為此標(biāo)準(zhǔn)形式 在MATLAB5 x以上版中 線性規(guī)劃問(wèn)題LinearProgramming已用函數(shù)linprog取代了MATLAB5 x版中的lp函數(shù) 當(dāng)然 由于版本的向下兼容性 一般說(shuō)來(lái) 低版本中的函數(shù)在7 0版中仍可使用 5 1線性優(yōu)化 函數(shù)linprog格式x linprog f A b 求minf xsub to線性規(guī)劃的最優(yōu)解 x linprog f A b Aeq beq 等式約束 若沒(méi)有不等式約束 則A b x linprog f A b Aeq beq lb ub 指定x的范圍 若沒(méi)有等式約束 則Aeq beq x linprog f A b Aeq beq lb ub x0 設(shè)置初值x0 x linprog f A b Aeq beq lb ub x0 options options為指定的優(yōu)化參數(shù) x fval linprog 返回目標(biāo)函數(shù)最優(yōu)值 即fval f x x lambda exitflag linprog lambda為解x的Lagrange乘子 x lambda fval exitflag linprog exitflag為終止迭代的錯(cuò)誤條件 x fval lambda exitflag output linprog output為關(guān)于優(yōu)化的一些信息 說(shuō)明 若exitflag 0表示函數(shù)收斂于解x exitflag 0表示超過(guò)函數(shù)估值或迭代的最大次數(shù) exitflag 0表示函數(shù)不收斂于解x 若lambda lower表示下界lb lambda upper表示上界ub lambda ineqlin表示不等式約束 lambda eqlin表示等式約束 lambda中的非0元素表示對(duì)應(yīng)的約束是有效約束 output iterations表示迭代次數(shù) output algorithm表示使用的運(yùn)算規(guī)則 output cgiterations表示PCG迭代次數(shù) 例5 1求下面的優(yōu)化問(wèn)題 min sub to 解 f 5 4 6 寫成行向量亦可 A 1 11 324 320 b 20 42 30 lb zeros 3 1 x fval exitflag output lambda linprog f A b lb 結(jié)果為 x 最優(yōu)解0 000015 00003 0000fval 最優(yōu)值 78 0000exitflag 收斂1 output iterations 6 迭代次數(shù)cgiterations 0algorithm lipsol 所使用規(guī)則lambda ineqlin 3x1double eqlin 0 x1double upper 3x1double lower 3x1double lambda ineqlinans 0 00001 50000 5000 lambda lower ans 1 00000 00000 0000表明 不等約束條件2和3以及第1個(gè)下界是有效的 請(qǐng)寫出下面線性規(guī)劃的Matlab程序 c 0 4 0 28 0 32 0 72 0 64 0 6 A 0 010 010 010 030 030 03 0 02000 0500 00 02000 050 000 03000 08 b 850 700 100 900 Aeq beq vlb 0 0 0 0 0 0 vub x fval linprog c A b Aeq beq vlb vub MATLAB求解優(yōu)化問(wèn)題的主要函數(shù) 優(yōu)化函數(shù)的輸入變量 優(yōu)化函數(shù)的輸出變量 5 2非線性優(yōu)化 5 2 1有約束的一元函數(shù)的最小值單變量函數(shù)求最小值的標(biāo)準(zhǔn)形式為sub to函數(shù)fminbnd格式x fminbnd fun x1 x2 返回自變量x在區(qū)間上函數(shù)fun取最小值時(shí)x值 fun為目標(biāo)函數(shù)的表達(dá)式字符串或MATLAB自定義函數(shù)的函數(shù)柄 x fminbnd fun x1 x2 options x fval fminbnd x fval exitflag fminbnd x fval exitflag output fminbnd 例5 2計(jì)算下面函數(shù)在區(qū)間 0 1 內(nèi)的最小值 解 x fval exitflag output fminbnd x 3 cos x x log x exp x 0 1 x 0 5223fval 0 3974exitflag 1 output iterations 9funcCount 9algorithm goldensectionsearch parabolicinterpolation 例5 3在 0 5 上求下面函數(shù)的最小值解 先自定義函數(shù) 在MATLAB編輯器中建立M文件為 functionf myfun x f x 3 2 1 保存為myfun m 然后在命令窗口鍵入命令 x fminbnd myfun 0 5 則結(jié)果顯示為 x 3 5 2 2無(wú)約束多元函數(shù)最小值多元函數(shù)最小值的標(biāo)準(zhǔn)形式為其中 x為向量 命令利用函數(shù)fminsearch求無(wú)約束多元函數(shù)最小值 函數(shù)fminsearch格式x fminsearch fun x0 x0為初始點(diǎn) fun為目標(biāo)函數(shù)的表達(dá)式字符串或MATLAB自定義函數(shù)的函數(shù)柄 x fminsearch fun x0 options options查optimset x fval fminsearch 最優(yōu)點(diǎn)的函數(shù)值 x fval exitflag fminsearch exitflag與單變量情形一致 x fval exitflag output fminsearch output與單變量情形一致 例5 4求的最小值點(diǎn) 解 X fminsearch 2 x 1 3 4 x 1 x 2 3 10 x 1 x 2 x 2 2 0 0 結(jié)果為X 1 00160 8335或在MATLAB編輯器中建立函數(shù)文件 functionf myfun x f 2 x 1 3 4 x 1 x 2 3 10 x 1 x 2 x 2 2 保存為myfun m 在命令窗口鍵入 X fminsearch myfun 0 0 或 X fminsearch myfun 0 0 結(jié)果為 X 1 00160 8335 5 2 3有約束的多元函數(shù)最小值非線性有約束的多元函數(shù)的標(biāo)準(zhǔn)形式為 sub to其中 x b beq lb ub是向量 A Aeq為矩陣 C x Ceq x 是返回向量的函數(shù) f x 為目標(biāo)函數(shù) f x C x Ceq x 可以是非線性函數(shù) 在MATLAB5 x中 它的求解由函數(shù)constr實(shí)現(xiàn) 函數(shù)fmincon格式x fmincon fun x0 A b x fmincon fun x0 A b Aeq beq x fmincon fun x0 A b Aeq beq lb ub x fmincon fun x0 A b Aeq beq lb ub nonlcon x fmincon fun x0 A b Aeq beq lb ub nonlcon options x fval fmincon x fval exitflag fmincon x fval exitflag output fmincon x fval exitflag output lambda fmincon x fval exitflag output lambda grad fmincon x fval exitflag output lambda grad hessian fmincon 參數(shù)說(shuō)明 fun為目標(biāo)函數(shù) 它可用前面的方法定義 nonlcon的作用是通過(guò)接受的向量x來(lái)計(jì)算非線性不等約束和等式約束分別在x處的估計(jì)C和Ceq 通過(guò)指定函數(shù)柄來(lái)使用 如 x fmincon myfun x0 A b Aeq beq lb ub mycon 先建立非線性約束函數(shù) 并保存為mycon m function C Ceq mycon x C 計(jì)算x處的非線性不等約束的函數(shù)值 Ceq 計(jì)算x處的非線性等式約束的函數(shù)值 lambda是Lagrange乘子 它體現(xiàn)哪一個(gè)約束有效 output輸出優(yōu)化信息 grad表示目標(biāo)函數(shù)在x處的梯度 hessian表示目標(biāo)函數(shù)在x處的Hessian值 例5 5求下面問(wèn)題在初始點(diǎn) 0 1 處的最優(yōu)解 min sub to 解 約束條件的標(biāo)準(zhǔn)形式為 sub to先在MATLAB編輯器中建立非線性約束函數(shù)文件 function c ceq mycon x c x 1 1 2 x 2 ceq 無(wú)等式約束 然后 在命令窗口鍵入如下命令或建立M文件 fun x 1 2 x 2 2 x 1 x 2 2 x 1 5 x 2 目標(biāo)函數(shù) x0 01 A 23 線性不等式約束 b 6 Aeq 無(wú)線性等式約束 beq lb x沒(méi)有下 上界 ub x fval exitflag output lambda grad hessian fmincon fun x0 A b Aeq beq lb ub mycon 則結(jié)果為x 34fval 13exitflag 1 解收斂 output iterations 2funcCount 9stepsize 1algorithm medium scale SQP Quasi Newton line search firstorderopt cgiterations lambda lower 2x1double x下界有效情況 通過(guò)lambda lower可查看 upper 2x1double x上界有效情況 為0表示約束無(wú)效 eqlin 0 x1double 線性等式約束有效情況 不為0表示約束有效 eqnonlin 0 x1double 非線性等式約束有效情況 ineqlin 2 5081e 008 線性不等式約束有效情況 neqnonlin 6 1938e 008 非線性不等式約束有效情況 grad 目標(biāo)函數(shù)在最小值點(diǎn)的梯度 1 0e 006 0 1776hessian 目標(biāo)函數(shù)在最小值點(diǎn)的Hessian值 1 0000 0 0000 0 00001 0000 5 2 4二次規(guī)劃問(wèn)題二次規(guī)劃問(wèn)題 quadraticprogramming 的標(biāo)準(zhǔn)形式為 sub to其中 H A Aeq為矩陣 f b beq lb ub x為向量其它形式的二次規(guī)劃問(wèn)題都可轉(zhuǎn)化為標(biāo)準(zhǔn)形式 MATLAB5 x版中的qp函數(shù)已被6 0版中的函數(shù)quadprog取代 函數(shù)quadprog格式x quadprog H f A b 其中H f A b為標(biāo)準(zhǔn)形中的參數(shù) x為目標(biāo)函數(shù)的最小值 x quadprog H f A b Aeq beq Aeq beq滿足等約束條件 x quadprog H f A b Aeq beq lb ub lb ub分別為解x的下界與上界 x quadprog H f A b Aeq beq lb ub x0 x0為設(shè)置的初值x quadprog H f A b Aeq beq lb ub x0 options options為指定的優(yōu)化參數(shù) x fval quadprog fval為目標(biāo)函數(shù)最優(yōu)值 x fval exitflag quadprog exitflag與線性規(guī)劃中參數(shù)意義相同 x fval exitflag output quadprog output與線性規(guī)劃中參數(shù)意義相同 x fval exitflag output lambda quadprog lambda與線性規(guī)劃中參數(shù)意義相同 例5 6求二次規(guī)劃的最優(yōu)解maxf x1 x2 x1x2 3sub tox1 x2 2 0解 化成標(biāo)準(zhǔn)形式 sub tox1 x2 2 在Matlab中實(shí)現(xiàn)如下 H 0 1 1 0 f 0 0 Aeq 11 b 2 x fval exitflag output lambda quadprog H f Aeq b 結(jié)果為 x 1 00001 0000 fval 1 0000exitflag 4output iterations 1algorithm large scale projectivepreconditionedconjugategradients firstorderopt 0cgiterations 1message Optimizationterminated localminimumfound thesolutionissingular lambda eqlin 1 0000ineqlin lower upper 5 3極小化極大 Minmax 問(wèn)題 sub to 其中 x b beq lb ub是向量 A Aeq為矩陣 C x Ceq x 和F x 是返回向量的函數(shù) F x C x Ceq x 可以是非線性函數(shù) 函數(shù)fminimax格式x fminimax fun x0 x fminimax fun x0 A b x fminimax fun x0 A b Aeq beq x fminimax fun x0 A b Aeq beq lb ub x fminimax fun x0 A b Aeq beq lb ub nonlcon x fminimax fun x0 A b Aeq beq lb ub nonlcon options x fval maxfval fminimax x fval maxfval exitflag fminimax x fval maxfval exitflag output fminimax x fval maxfval exitflag output lambda fminimax 例5 7求下列函數(shù)最大值的最小化問(wèn)題其中 解 先建立目標(biāo)函數(shù)文件 并保存為myfun m functionf myfun x f 1 2 x 1 2 x 2 2 48 x 1 40 x 2 304 f 2 x 1 2 3 x 2 2 f 3 x 1 3 x 2 18 f 4 x 1 x 2 f 5 x 1 x 2 8 然后 在命令窗口鍵入命令 x0 0 1 0 1 初始值 x fval fminimax myfun x0 結(jié)果為 x 4 00004 0000fval 0 0000 64 0000 2 0000 8 0000 0 0000 5 4曲線擬合與插值 在大量的應(yīng)用領(lǐng)域中 人們經(jīng)常面臨用一個(gè)解析函數(shù)描述數(shù)據(jù) 通常是測(cè)量值 的任務(wù) 對(duì)這個(gè)問(wèn)題有兩種方法 插值 在插值法里 數(shù)據(jù)假定是正確的 要求以某種方法描述數(shù)據(jù)點(diǎn)之間所發(fā)生的情況 曲線擬合 曲線擬合或回歸是人們?cè)O(shè)法找出某條光滑曲線 它最佳地?cái)M合數(shù)據(jù) 但不必要經(jīng)過(guò)任何數(shù)據(jù)點(diǎn) 標(biāo)有 o 的是數(shù)據(jù)點(diǎn) 連接數(shù)據(jù)點(diǎn)的實(shí)線描繪了線性內(nèi)插 虛線是數(shù)據(jù)的最佳擬合 曲線擬合的兩個(gè)基本問(wèn)題 1 最佳擬合意味著什么 2 應(yīng)該用什么樣的曲線 可用許多不同的方法定義最佳擬合 并存在無(wú)窮數(shù)目的曲線 當(dāng)最佳擬合被解釋為在數(shù)據(jù)點(diǎn)的最小誤差平方和 且所用的曲線限定為多項(xiàng)式時(shí) 那么曲線擬合是相當(dāng)簡(jiǎn)捷的 數(shù)學(xué)上 稱為多項(xiàng)式的最小二乘曲線擬合 1Matlab曲線擬合 虛線和標(biāo)志的數(shù)據(jù)點(diǎn)之間的垂直距離是在該點(diǎn)的誤差 對(duì)各數(shù)據(jù)點(diǎn)距離求平方 并把平方距離全加起來(lái) 就是誤差平方和 這條虛線是使誤差平方和盡可能小的曲線 即是最佳擬合 最小二乘這個(gè)術(shù)語(yǔ)僅僅是使誤差平方和最小的省略說(shuō)法 Matlab曲線擬合和插值命令 x 0 1 2 3 4 5 6 7 8 91 y 4471 9783 286 167 087 347 669 569 489 3011 2 n 2 p polyfit x y n ezplot 9 8108 x x 20 1293 x 0 0317 二次多項(xiàng)式系數(shù) 既是p的輸出 該命令為畫出擬合曲線 xi linspace 0 1 100 x axisdataforplottingz polyval p xi 求出多項(xiàng)式在xi點(diǎn)處的取值 plot x y o x y xi z xlabel x ylabel y f x title SecondOrderCurveFitting 多項(xiàng)式階次的選擇是任意的 兩點(diǎn)決定一直線或一階多項(xiàng)式 三點(diǎn)決定一個(gè)平方或2階多項(xiàng)式 按此進(jìn)行 n 1數(shù)據(jù)點(diǎn)唯一地確定n階多項(xiàng)式 于是 在上面的情況下 有11個(gè)數(shù)據(jù)點(diǎn) 我們可選一個(gè)高達(dá)10階的多項(xiàng)式 然而 高階多項(xiàng)式給出很差的數(shù)值特性 人們不應(yīng)選擇比所需的階次高的多項(xiàng)式 此外 隨著多項(xiàng)式階次的提高 近似變得不夠光滑 因?yàn)檩^高階次多項(xiàng)式在變零前 可多次求導(dǎo) 原始數(shù)據(jù)標(biāo)以 o 2階曲線擬合是虛線 10階擬合是實(shí)線 注意 在10階擬合中 在左邊和右邊的極值處 數(shù)據(jù)點(diǎn)之間出現(xiàn)大的紋波 當(dāng)企圖進(jìn)行高階曲線擬合時(shí) 這種紋波現(xiàn)象經(jīng)常發(fā)生 顯然 越多就越好 的觀念在這里不適用 2插值命令 插值定義為對(duì)數(shù)據(jù)點(diǎn)之間函數(shù)的估值方法 這些數(shù)據(jù)點(diǎn)是由某些集合給定 當(dāng)人們不能很快地求出所需中間點(diǎn)的函數(shù)值時(shí) 插值是一個(gè)有價(jià)值的工具 例如 當(dāng)數(shù)據(jù)點(diǎn)是某些實(shí)驗(yàn)測(cè)量的結(jié)果或是過(guò)長(zhǎng)的計(jì)算過(guò)程時(shí) 就有這種情況 最簡(jiǎn)單插值的例子是MATLAB的作圖 按缺省 MATLAB用直線連接所用的數(shù)據(jù)點(diǎn)以作圖 這個(gè)線性插值猜測(cè)中間值落在數(shù)據(jù)點(diǎn)之間的直線上 當(dāng)然 當(dāng)數(shù)據(jù)點(diǎn)個(gè)數(shù)的增加和它們之間距離的減小時(shí) 線性插值就更精確 例如 x1 linspace 0 2 pi 60 x2 linspace 0 2 pi 6 plot x1 sin x1 x2 sin x2 xlabel x ylabel sin x title LinearInterpolation 一個(gè)在數(shù)據(jù)點(diǎn)之間用60個(gè)點(diǎn) 它比另一個(gè)只用6個(gè)點(diǎn)更光滑和更精確 如曲線擬合一樣 插值要作決策 根據(jù)所作的假設(shè) 有多種插值 而且 可以在一維以上空間中進(jìn)行插值 即如果有反映兩個(gè)變量函數(shù)的插值 z f x y 那么就可在x之間和在y之間 找出z的中間值進(jìn)行插值 MATLAB在一維函數(shù)interp1和在二維函數(shù)interp2中 提供了許多的插值選擇 為了說(shuō)明一維插值 考慮下列問(wèn)題 12小時(shí)內(nèi) 一小時(shí)測(cè)量一次室外溫度 數(shù)據(jù)存儲(chǔ)在兩個(gè)MATLAB變量中 hours 1 12 indexforhourdatawasrecordedtemps 589152529313022252724 recordedtemperaturesplot hours temps hours temps viewtemperaturestitle Temperature xlabel Hour ylabel DegreesCelsius 正如上圖看到的 MATLAB畫出了數(shù)據(jù)點(diǎn)線性插值的直線 為了計(jì)算在任意給定時(shí)間的溫度 人們可試著對(duì)可視的圖作解釋 另外一種方法 可用函數(shù)interp1 t interp1 hours temps 9 3 estimatetemperatureathour 9 3t 22 9000t interp1 hours temps 4 7 estimatetemperatureathour 4 7t 22t interp1 hours temps 3 26 57 111 7 findtempatmanypoints t 10 200030 000030 900024 9000 interp1的缺省用法是由interp1 x y xo 來(lái)描述 這里x是獨(dú)立變量 橫坐標(biāo) y是應(yīng)變量 縱坐標(biāo) xo是進(jìn)行插值的一個(gè)數(shù)值數(shù)組 另外 該缺省的使用假定為線性插值 若不采用直線連接數(shù)據(jù)點(diǎn) 我們可采用某些更光滑的曲線來(lái)擬合數(shù)據(jù)點(diǎn) 最常用的方法是用一個(gè)3階多項(xiàng)式 即3次多項(xiàng)式 來(lái)對(duì)相繼數(shù)據(jù)點(diǎn)之間的各段建模 每個(gè)3次多項(xiàng)式的頭兩個(gè)導(dǎo)數(shù)與該數(shù)據(jù)點(diǎn)相一致 這種類型的插值被稱為3次樣條或簡(jiǎn)稱為樣條 函數(shù)interp1也能執(zhí)行3次樣條插值 t interp1 hours temps 9 3 spline estimatetemperatureathour 9 3t 21 8577t interp1 hours temps 4 7 spline estimatetemperatureathour 4 7t 22 3143t interp1 hours temps 3 26 57 111 7 spline t 9 673430 042731 175525 3820 interp1二個(gè)強(qiáng)約束 1 人們不能要求有獨(dú)立變量范圍以外的結(jié)果 例如 interp1 hours temps 13 5 導(dǎo)致一個(gè)錯(cuò)誤 因?yàn)閔ours在1到12之間變化 2 獨(dú)立變量必須是單調(diào)的 即獨(dú)立變量在值上必須總是增加的或總是減小的 二維插值是基于與一維插值同樣的基本思想 然而 正如名字所隱含的 二維插值是對(duì)兩變量的函數(shù)z f x y 進(jìn)行插值 3 非線性數(shù)據(jù) 曲線 擬合非線性曲線擬合是已知輸入向量xdata和輸出向量ydata 并且知道輸入與輸出的函數(shù)關(guān)系為ydata F x xdata 但不知道系數(shù)向量x 今進(jìn)行曲線擬合 求x使得下式成立 函數(shù)lsqcurvefit格式x lsqcurvefit fun x0 xdata ydata x lsqcurvefit fun x0 xdata ydata lb ub options x resnorm lsqcurvefit x resnorm residual exitflag output lambda lsqcurvefit x resnorm residual exitflag output lambda jacobian lsqcurvefit resnorm sum fun x xdata ydata 2 即在x處殘差的平方和 residual fun x xdata ydata 即在x處的殘差 exitflag為終止迭代的條件 output為輸出的優(yōu)化信息 lambda為解x處的Lagrange乘子 jacobian為解x處擬合函數(shù)fun的jacobian矩陣 例5 8求解如下最小二乘非線性擬合問(wèn)題已知輸入向量xdata和輸出向量ydata 且長(zhǎng)度都是n 擬合函數(shù)為即目標(biāo)函數(shù)為其中 初始解向量為x0 0 3 0 4 0 1 解 先建立擬合函數(shù)文件 并保存為myfun mfunctionF myfun x xdata F x 1 xdata 2 x 2 sin xdata x 3 xdata 3 然后給出數(shù)據(jù)xdata和ydata xdata 3 67 79 34 18 62 81 37 910 05 4 ydata 16 5150 6263 124 7208 59 92 7163 9325 054 3 x0 10 10 10 初始估計(jì)值 x resnorm lsqcurvefit myfun x0 xdata ydata 結(jié)果為 Optimizationterminatedsuccessfully RelativefunctionvaluechangingbylessthanOPTIONS TolFunx 0 22690 33850 3021resnorm 6 2950 4非線性最小二乘非線性最小二乘 非線性數(shù)據(jù)擬合 的標(biāo)準(zhǔn)形式為 其中 L為常數(shù)設(shè)則目標(biāo)函數(shù)可表達(dá)為其中 x為向量 F x 為函數(shù)向量 函數(shù)lsqnonlin格式x lsqnonlin fun x0 x0為初始解向量 fun為 i 1 2 m fun返回向量值F 而不是平方和值 平方和隱含在算法中 fun的定義與前面相同 x lsqnonlin fun x0 lb ub options options為指定優(yōu)化參數(shù) 若x沒(méi)有界 則lb ub x resnorm lsqnonlin resnorm sum fun x 2 即解x處目標(biāo)函數(shù)值 x resnorm residual lsqnonlin residual fun x 即解x處fun的值 例5 9求下面非線性最小二乘問(wèn)題初始解向量為x0 0 3 0 4 解 先建立函數(shù)文件 并保存為myfun m 由于lsqnonlin中的fun為向量形式而不是平方和形式 因此 myfun函數(shù)應(yīng)由建立 k 1 2 10functionF myfun x k 1 10 F 2 2 k exp k x 1 exp k x 2 然后調(diào)用優(yōu)化程序 x0 0 30 4 x resnorm lsqnonlin myfun x0 x 0 25780 2578resnorm 求目標(biāo)函數(shù)值124 3622 實(shí)際問(wèn)題求解 任務(wù)分配問(wèn)題 某車間有甲 乙兩臺(tái)機(jī)床 可用于加工三種工件 假定這兩臺(tái)車床的可用臺(tái)時(shí)數(shù)分別為800和900 三種工件的數(shù)量分別為400 600和500 且已知用三種不同車床加工單位數(shù)量不同工件所需的臺(tái)時(shí)數(shù)和加工費(fèi)用如下表 問(wèn)怎樣分配車床的加工任務(wù) 才能既滿足加工工件的要求 又使加工費(fèi)用最低 解設(shè)在甲車床上加工工件1 2 3的數(shù)量分別為x1 x2 x3 在乙車床上加工工件1 2 3的數(shù)量分別為x4 x5 x6 可建立以下線性規(guī)劃模型 f 1391011128 A 0 41 110000000 51 21 3 b 800 900 Aeq 100100010010001001 beq 400600500 vlb zeros 6 1 vub x fval linprog f A b Aeq beq vlb vub 程序詳解 見example5 m Theend