用有限的資源得到最大的收益,或是用最少的資源得到足夠的收益,是作業研究(Operations Research)領域當中最基礎的兩個議題,也是高中或大學時期所學到的線性規劃(Linear Programming)一定會接觸到的命題。就算是在日常生活也很有可能遇到這類問題,一般而言我們遇到的問題也不會太複雜,都是線性規劃就可以解決的問題。這篇文章將以筆者自身在遊玩蔚藍檔案日版時,對於復刻情人節活動的資源規劃過程為例,示範線性規劃的問題整理方法,以及如何使用電腦程式的輔助算出解答。
前言
一切都是為了與所有人約會,以及合歓垣フブキ!

因為筆者是2022年8月才開始遊玩蔚藍檔案,沒有參與過情人節活動,正好遇到今年1月的復刻,也不知道下次再復刻是何時,所以想要一口氣將目前手上的角色都約會過一次。加上筆者實際把フブキ拿到手之後,發覺我好像滿喜歡她這種感覺的,想說反正都要花大筆成本拿所有角色的巧克力了,就順便定下目標:將フブキ升到5星。除此之外,我還要把商店的物品全部搬空。
可想而知要達到這些目標,我需要每天花非常多AP。加上關卡獎勵的設計,我沒有辦法簡單地得到「全部刷這關就好!」的結論。這次活動還有個讓事情變得更複雜的系統:每日花費一定AP之後會觸發狐坂ワカモ的突襲事件,而這個突襲事件的獎勵特別豐厚,規劃問題時不能將這點排除在外。面對這種複雜的命題,正是線性規劃的用武之地。
不過要注意的是,上頭提到的獎勵都是確定性的,關卡掉落的通貨是確定數量的,突襲事件的掉落也是固定的。線性規劃沒辦法幫你規劃隨機性掉落的物品,比如オーパーツ,如果有這方面需求的話那先刷到滿足的數量再來做規劃會比較好。
問題架構
問題敘述
活動日期
在問題思考的當下,我總共還有12個完整的遊戲天+最後1天維修前,總共13天的時間。
最後1天維修前的時間雖然比較短,但以一天回復6次AP共720AP+戰術對抗戰商店4次90AP共360AP,加總共1080AP來看,已經足夠觸發當日所有突襲事件,所以大可以把最後1天當成完整的一天來看。
同時我們也可以規定一天最多能花費的AP數量,避免算出什麼第一天花9000AP之後每天不用花的離譜結果。
關卡範圍
活動關卡有12關,其中9~12關是難度最高但效率也最高的,所以只考慮這四關的AP消耗與掉落通貨。
這些關卡會掉落以下四種通貨:PT(活動分數)、問題児の反省文、押収された違法駄菓子、環境美化活動用ゴミ袋。
通貨定義
PT作為通貨的功能是換取約會邀請,550PT換一張,所以需求量即為550乘以學生數量。除此之外累積PT點也能拿一些獎勵,累積10000點即可拿完。由於我有許多角色等著約,累積的PT點必然是會超過10000點的,就算計算問題的當下還沒拿到這些獎勵,還是得提前把這些獎勵納入考量(參考)。
其它三種通貨可以在商店裡頭換取獎勵。我們可以簡單地算出我們需要多少通貨才能把商店中有限次交易的物品搬空。同時也存在通貨間的無限次交換,包含以5個反省文換違法糖果、以5個違法糖果換垃圾袋,以及以200個垃圾袋換フブキ的神名文字(將神名文字也視為通貨的原因請參照下一段落)。
這邊要注意的是,計算帳面需求量時不要將通貨之間的交換納入考量。比如我今天規定要「100個反省文與100個違法糖果」,我不能擅自把規定改寫成「600個反省文(其中500個換成100個違法糖果)」,因為這等同規定了「100個違法糖果裡頭有多少是從反省文換來的」這件事。但這根本不是我們在意的事情,我們在設計問題時只在乎「確實拿到100個反省文用來搬空商店」與「確實拿到100個違法糖果用來搬空商店」這種帳面需求量。至於會有多少額外的反省文讓我們換成違法糖果,那是線性規劃該去處理的事情。
神名文字
フブキ的神名文字雖然不是上述列出的通貨,本質上神名文字也確實不是通貨(不能拿來換其它物品),但在計算上把神名文字當成一種通貨是有好處的。因為其它四種通貨我們都已經有其帳面需求量,而神名文字也有自己的需求量(比如我要從3星升到5星需要220顆)、在關卡的掉落量是固定的(一般關卡不掉、突襲關卡固定1顆),而且也有將垃圾袋換成神名文字的無限次交換管道。為了描述方便,我們也將神名文字也視為一種通貨。
這邊也有需要注意的點,在計算神名文字的需求量時務必將「累積PT」得到的神名文字獎勵給考慮進去。例如在規劃問題的當下還需要220顆神名文字,但累積PT獎勵還有20顆沒有拿到,那麼需求量就是200顆,因為在刷關的過程中,必然會把那20顆獎勵也一併拿到手。除非你的需求真的很低,導致總PT打不到6000……這樣的話也沒必要用線性規劃徒增麻煩吧。
突襲事件
每天——包含最後1天換日後到維修前的期間——在活動花費「20、50、100、180、300、540、780、1020」AP時各會得到一張「ワカモの逮捕状」,並且觸發ワカモ突襲事件。打一次可以獲得固定數量的獎勵,包含上述所有種類的通貨。
如果是像我這樣今年才第一次參加情人節活動的新老師,為了與所有學生約會需要付出比較大的代價,那麼每日1020AP基本上是必要的,可以把每日8次的突襲事件當成前提,總獎勵量直接折抵上面各個通貨的需求量,以簡化問題。不過這篇文章筆者決定依然將每日突襲事件的發生次數納入考量,除了保留彈性讓去年參加過的老手玩家也能夠適當規劃資源之外,也能示範線性規劃如何解決這種進階命題。
目標
我們的目標是最小化這13天總共花費在活動上面的AP量。
至此,我們已經將問題架構用中文稍微爬梳過一次了。接下來我們將會把這些敘述抽象化並轉成數學式。
抽象化與數學式
參數與變數
首先要解釋幾個觀念:參數(parameters)與變數(variables)。
「參數」是在計算過程中會引用到的數值,這些值不會因為計算而改變,一般人可以簡單把它們理解成是「固定的值」或。在這個題目裡頭有這些值:總共天數(13)、關卡9掉落的PT量、關卡9消耗的AP(20)、花費多少垃圾袋交換1顆神名文字(200)等等。
關於「參數」與「常數」的差異……
嚴格來說裡頭有些數值不是我們能調整的,比如關卡9消耗的AP量以及交換匯率等;其它參數的值則是會視自身情況調整,比如總共天數以及關卡掉落的PT量等。前者有時我們會稱之為「常數」來與「參數」做區隔,這種值在列式時可以直接寫出數字而不用寫成代號。但為了列式子方便,有些「常數」亦會寫成代號的形式,比如「關卡9花費的AP」與其寫死成$20$,寫成$Apcost_9$更方便用於$\sum_s Apcost_s$之類的情境。在之後的文章不會再強調「參數」與「常數」的差異性。
「變數」則是在計算過程中會變化的值,也是我們想要「算出來」的值。在這個題目裡頭大概就是:今天這個關卡要「刷幾次」這類的。
線性規劃之所以「線性」,就是因為我們的數學式都是各個變數的線性組合,不會有變數乘以變數(二次方)或者對變數開平方根的運算。大部分生活中的情境都能用線性組合來表述,這邊的題目也不例外。
集合
為了表達「關卡9花費的AP」,我們可能會寫出$Ap\_cost_9$這個代號,下標的9換成其它數值即可代換成其它關卡的AP花費量。我們只考慮關卡9~12,所以我們會有$Ap\_cost_9$、$Ap\_cost_{10}$、$Ap\_cost_{11}$與$Ap\_cost_{12}$這四個代號。為了更方便整理這些代號,我們會定義一個「集合」$S := \{9, 10, 11, 12\}$,而這些代號就會變成$Ap\_cost_s, s \in S$。
一開始可能還不知道要定義哪些集合才夠。一般而言在抽象化的過程裡,定義參數與變數時才會找出需要先定義哪些集合。但因為這又必須定義在引用到它們的參數或變數之前,所以筆者先把結論列在最前頭:
| 名稱 | 值 | 意義 |
| $D$ | $1..Days$ | 第1天到第$Days$天(Days為參數) |
| $M$ | $\{20, 50, 100, 180, 300, 540, 780, 1020\}$ | 觸發突襲的AP量(M取自milestone) |
| $S$ | $9..12$ | 考量的關卡(S取自stage) |
| $C$ | $\{point, paper, candy, bag, moji\}$ | 通貨種類(C取自current) 分別代表PT、反省文、違法的糖果、垃圾袋與神名文字 |
參數
我們慢慢以先前問題敘述的段落順序來列出參數,參數的代號一律以大寫開頭表示。參數的值不需要急著指定。下列出現的下標小寫字母分別對應到大寫字母的集合,比如$Ap\_cost_s$的$s$即為集合$S$的元素。
活動日期
活動天數定義成$Days$,每日最多花費的AP量定義成$Max\_ap$。
關卡範圍
關卡$s$花費的AP定義成$Ap\_cost_s$,關卡$s$掉落的通貨$c$數量為$Stage\_get_{s, c}$。
通貨定義+神名文字
各種通貨$c$的帳面需求量為$Need_c$。
對於通貨間交易,反省文與違法糖果的交換匯率為$Rate\_paper2candy$、違法糖果與垃圾袋的交換匯率為$Rate\_candy2bag$,垃圾袋與神名文字的交換匯率則為$Rate\_bag2moji$。我們不用$Rate_{c, c}$來表示是因為總共25種可能性當中我們只會用到3種,沒有必要徒增麻煩。
突襲事件
刷一次突襲關卡可以拿到的通貨$c$為$M\_get_c$。
統整下來我們有下列參數,值為空的部分代表需要依照自身情況調整:
| 名稱 | 值 | 意義 |
| $Days$ | 總天數 | |
| $Max\_ap$ | 每日花費的AP量限制 | |
| $Ap\_cost_s$ | $9: 20$ $10: 20$ $11: 20$ $12: 20$ | 關卡$s$花費的AP |
| $Stage\_get_{s, c}$ | 關卡$s$掉落的通貨$c$數量 | |
| $Need_c$ | 通貨$c$的帳面需求量 | |
| $Rate\_paper2candy$ | $5$ | 5個反省文換1個違法糖果 |
| $Rate\_candy2bag$ | $5$ | 5個違法糖果換1個垃圾袋 |
| $Rate\_bag2moji$ | $200$ | 200個垃圾袋換1顆神名文字 |
| $M\_get_c$ | $point: 120$ $paper: 40$ $candy: 35$ $bag: 30$ $moji: 1$ | 突襲關卡的獎勵 |
變數與條件式
我們的問題有兩大因素要考量,一個是我們的「目標」最小化AP總花費,一個是「條件」所有通貨都獲得至少一定的數量去換取約會機會/搬空商店/升フブキ的星級。變數的代號都是以小寫字母為開頭。
本題出現的變數都是非負整數,也就是大於等於$0$的整數。這件事將會事先告訴我們的解題器。
目標
我們的目標是最小化AP總花費。因為後續還要考慮每天觸發突襲的有無,所以將「AP總花費」從日期的角度切入,即為「每天各自花費的AP總量」的加總,定義「第$d$天花費的AP總量」為$ap_d$。於是我們的目標式為:
\[ \textbf{minimize} \sum_{d \in D} ap_d \]
我們還可以細分第幾天花費在某一關多少AP,但我們直接定義「第$d$天打關卡$s$的次數」為$time_{d, s}$即可,兩者只差在乘以$Ap\_cost_s$。於是可以得到每天花費的AP總量與$time_{d, s}$的關係式:
\[ \forall d \in D. ap_d = \sum_{s \in S} Ap\_cost_s \times time_{d, s} \]
然後不要忘記,我們還要限制每天最多花多少AP:
\[ \forall d \in D. ap_d \le Max\_ap \]
通貨需求
其實有了$time_{d, s}$之後,我們就可以算出打關卡直接得到的通貨總量。不過要得出與需求量$Need_c$的關係式,我們還需要考慮通貨的另外兩種入手手段:一種是ワカモ突襲關卡,一種是通貨間的交換。
為了得出突襲關卡獲得的通貨總量,我們需要得知突襲關卡被觸發了多少次。我們可以定義每天觸發幾次突襲關卡,不過這邊我們直接定義「第$d$天花費達到$m$ AP與否」為$to\_ap_{d, m}$,它的值是$0$或$1$,達到的話值為$1$,即代表觸發1次。加總起來就可以得到第$d$天觸發的次數。
當觸發次數變多,很自然地得到的通貨量也會變多。但我們必須告訴線性規劃解題者說:觸發次數不見得越多越好,因為它需要付出AP作為代價,才能避免最後得到的解答將$to\_ap_{d, m}$全都調成$1$但是花費總AP對不上的狀況。「如果這個開關為1則某件事情要成立」是線性規劃常見的一種狀況,我們會用下列這種方法限制:
\[ \forall d \in D, m \in M. ap_d \ge m \times to\_ap_{d, m} \]
例如$to\_ap_{1, 540}=1$的話,則$ap_d$就要大於或等於$540$。
雖然為$0$的話$ap_d$也會限制成要大於等於$0$,但$ap_d$本來就是非負整數,所以相當於沒有限制。
有了打關卡的次數,我們就可以得到打關卡直接拿到的通貨量。因為我們還有通貨交換要考量,所以打關卡拿到的通貨量我們定為$base\_get_c$,透過通貨交易得到的通貨$c$的總量定義成$ext\_get_c$,兩者的加總為$total\_get_c$:
\[ \forall c \in C. base\_get_c + ext\_get_c = total\_get_c \]
這個加總當然必須至少達到我們對此種通貨的帳面需求量:
\[ \forall c \in C. total\_get_c \ge Need_c \]
總共入手量減去帳面需求量就是「多出來」的通貨,白話來說就是我已經把商店裡頭這通貨能換出來的物品搬空了,多出來的通貨才會——也只能——拿去換成其它種類的通貨。
先把目光轉回推關的收入。打關卡直接拿到的通貨量分成一般關卡與突襲關卡:
- 一般關卡:$\forall c \in C. \sum_{s \in S}(\sum_{d \in D} time_{d, s} \times Stage\_get_{s, c})$
每關卡刷的總次數乘以在此關卡得到通貨$c$的數量。 - 突襲關卡:$\forall c \in C. \sum_{d \in D} \sum_{m \in M} to\_ap_{d, m} \times M\_get_c$
觸發突襲的次數乘以突襲關卡得到通貨$c$的數量。
加起來就變成:
\[ \forall c \in C. base\_get_c = \sum_{s \in S}(\sum_{d \in D} time_{d, s} \times Stage\_get_{s, c}) + \sum_{d \in D} \sum_{m \in M} to\_ap_{d, m} \times M\_get_c \]
最後只剩下「通貨間交易」要解決了。首先我們已經知道有些通貨是沒有辦法被換出來的,包括PT與反省文,於是:
\[ ext\_get_{point} = 0 \]
\[ ext\_get_{paper} = 0 \]
也可以寫成$base\_get = total\_get$的形式,總之只要讓解題器知道這件事即可。
接下來我們先聚焦在「5反省文換1違法糖果」的交換途徑。
我們拿到的所有反省文是$total\_get_{paper}$,去掉帳面需求量$Need_{paper}$之後就是我們多出來的反省文。這些多出來的量:$total\_get_{paper} – Need_{paper}$才是我們能夠交換成違法糖果的。所以可以透過交易獲得的違法糖果數量是:
\[ ext\_get_{candy} = \floor{\frac{total\_get_{paper} – Need_{paper}}{Rate\_paper2candy}} \]
但是需要注意!$y = \floor{x}$並不是線性函數,所以不能用在「線性」規劃裡頭。那怎麼辦?
之所以會使用向下取整,就是因為$total\_get_{paper} – Need_{paper}$不見得是$5$(也就是$Rate\_paper2candy$)的倍數嘛。所以我們只需要有一個方法將分子變成$5$的倍數就好了……減去那些不滿$5$個的數量就可以了。這些「多出來」的反省文因為不足$5$個所以連換都沒得換,我們令這些沒有用的多餘反省文的數量叫做$dummy_{paper}$,這樣就可以改寫成:
\[ ext\_get_{candy} = \frac{total\_get_{paper} – Need_{paper} – dummy_{paper}}{Rate\_paper2candy} \]
因為解題器知道上面出現的這些變數都是非負整數,所以解題器不會亂歪一些奇怪的小數。
也別忘了限制$dummy_{paper}$的數值為$0$到$4$。因為已知$dummy_{paper}$是非負整數,我們只需要寫:
\[ dummy_{paper} \le Rate\_paper2candy\ -\ 1 \]
為何不寫「小於」
一般線性規劃的限制式是不允許寫$<$與$>$的。
假設有個變數$x$的範圍是$x < 5$,而目標是讓$x$最大化,那線性規劃就會一直想趨近$5$但又找不到終點。這是源於實數的稠密性。
你可能會問:啊現在這邊的例子左右兩邊都是整數,不會有上述的問題發生,為什麼還是不開放我們寫$<$呢?解題器在讀到這裡時,需要特地翻上下文才知道左右兩邊都是整數。如果只有這種等級的查詢還好,如果是類似費波那契數列一般項$\frac{\sqrt{5}}{5} \cdot [(\frac{1 + \sqrt{5}}{2})^2 – (\frac{1 – \sqrt{5}}{2})^2]$那樣,無法直接判斷是整數的式子怎麼辦?縱觀來說這種情況只是特例,而且我們列式時就可以用$\le$與$\ge$改寫成等價的式子,所以沒有特地允許$<$與$>$出現的必要。
可不可以直接說「$dummy_{paper}$是介於0到4之間的整數」而不是用限制式定義上限
可以。
在數學上這兩者其實沒差,因為只需要知道變數是整數,「非負整數」的「非負」部分就可以用$\ge 0$這樣的線性條件式限制出來。
但之所以會有這個疑慮,是因為在解題器運作層面上,「事先知道值域」與「看到條件式再去限縮」是有很大差異的,事先把值域定義好可以簡化一些運算……之後使用解題器時再說吧,總之這邊我們維持這樣的寫法即可。
另外兩個交易也能寫出來:
\[ ext\_get_{bag} = \frac{total\_get_{candy} – Need_{candy} – dummy_{candy}}{Rate\_candy2bag} \]
\[ dummy_{candy} \le Rate\_candy2bag\ -\ 1 \]
\[ ext\_get_{moji} = \frac{total\_get_{bag} – Need_{bag} – dummy_{bag}}{Rate\_bag2moji} \]
\[ dummy_{bag} \le Rate\_bag2moji\ -\ 1 \]
統整下來我們有下列變數:
| 名稱 | 型態/值域 | 意義 |
| $ap_d$ | 非負整數 | 第$d$天花費的AP總量 |
| $time_{d, s}$ | 非負整數 | 第$d$天打關卡$s$的次數 |
| $to\_ap_{d, m}$ | 0或1 | 第$d$天花費達到$m$ AP與否 |
| $base\_get_c$ | 非負整數 | 打關卡直接得到通貨$c$的總量 |
| $ext\_get_c$ | 非負整數 | 透過交易得到的通貨$c$的總量 |
| $total\_get_c$ | 非負整數 | 上兩者的加總 |
| $dummy_c$ | 非負整數 | 通貨交易的餘數 |
有關$dummy_c$
與$ext\_get_{point}$不同,$dummy_{point}$這類沒用到的變數不需要特別限制它等於$0$。因為在我們的列式當中,$ext\_get_{point}$會影響$total\_get_{point}$,如果沒有對$ext\_get_{point}$加以限制,解題器會直接令$ext\_get_{point}$為很大的數字以滿足$total\_get_{point}$,那樣算出來的東西就是錯的、沒意義的。但$dummy_{point}$則沒有參與任何限制式,所以它的值就算被亂標也不會影響。當然你大可以直接個別定義三個會用到的$dummy$,而不是像我這樣對所有通貨$c$都宣告一個。
以及這些式子:
| 式子 | 意義 |
| $\textbf{minimize} \sum_{d \in D} ap_d$ | 目標式:最小化AP總花費 |
| $\forall d \in D. ap_d = \sum_{s \in S} Ap\_cost_s \times time_{d, s}$ | 每日花費AP的總量 |
| $\forall d \in D. ap_d \le Max\_ap$ | 每天最多花多少AP |
| $\forall d \in D, m \in M. ap_d \ge m \times to\_ap_{d, m}$ | 花費一定AP才觸發突襲 |
| $\forall c \in C. base\_get_c + ext\_get_c = total\_get_c$ | 總共得到的通貨 |
| $\forall c \in C. total\_get_c \ge Need_c$ | 滿足需求量 |
| $\begin{split} \forall c \in C. base\_get_c = & \sum_{s \in S}(\sum_{d \in D} time_{d, s} \times Stage\_get_{s, c}) + \\ & \sum_{d \in D} \sum_{m \in M} to\_ap_{d, m} \times M\_get_c \end{split}$ | 從關卡獲得的通貨 |
| $ext\_get_{point} = 0$ $ext\_get_{paper} = 0$ | 無法透過交易獲得的通貨 |
| $ext\_get_{candy} = \frac{total\_get_{paper} – Need_{paper} – dummy_{paper}}{Rate\_paper2candy}$ $dummy_{paper} \le Rate\_paper2candy\ -\ 1$ | 反省文換違法糖果 |
| $ext\_get_{bag} = \frac{total\_get_{candy} – Need_{candy} – dummy_{candy}}{Rate\_candy2bag}$ $dummy_{candy} \le Rate\_candy2bag\ -\ 1$ | 違法糖果換垃圾袋 |
| $ext\_get_{moji} = \frac{total\_get_{bag} – Need_{bag} – dummy_{bag}}{Rate\_bag2moji}$ $dummy_{bag} \le Rate\_bag2moji\ -\ 1$ | 垃圾袋換神名文字 |
我們總算把問題整理完畢了,接下來就是告訴解題器這些東西,並且讓它去解題。
解題
我們要召喚我們的解題軟體了:SCIP(Solving Constraint Integer Programs),這是一款免費非商用的軟體,提供一個架構並依照不同的問題種類(包含基本的非線性規劃)調用內建的解題器進行解題。在這邊我們不需要知道太多內部細節,只需要用SCIP看得懂的方式告訴它題目即可。
寫出一個SCIP看得懂的題目有很多種方法,SCIP自己即支援多種不同的線性規劃或非線性規劃的撰寫格式,這些格式或有精細或有方便,適合用在不同的場合。這邊我們並不要求格式精簡,以「方便撰寫」為目標,我推薦使用ZIMPL(Zuse Institut Mathematical Programming Language)格式,寫起來的體驗最接近AMPL提供的格式。
另外一種方式則是使用其它API調用SCIP,比如Python的PySCIPOpt,它們會自動翻譯成SCIP讀得懂的格式並餵給SCIP,好處是可以使用自己熟悉的程式語言進行開發。當然讀者不見得會寫程式,所以不知道這是多棒的一件事。但這代表筆者可以輕鬆地用Flask架一個網站服務並呼叫SCIP解題,就可以在這篇文章底下弄一個表格讓讀者進行線上解題……
下面的題目撰寫將會用分頁的方式分成ZIMPL與PySCIPOpt兩種代碼進行示範。
輸入集合與參數
在ZIMPL裡頭,用set關鍵字去宣告集合、param宣告參數,後面加上各自的名稱,並且用:=去指定值。
集合值的語法是如同數學式那樣使用{}包圍,元素以逗號分隔,而連續整數可以用..標記,例如{1..3}等同{1, 2, 3}。
每個句子的尾巴要用分號做結尾。
像$point$這樣的字串需要用雙引號包圍。#則代表後面接的是註解。
param Days := 13; #因為D集合定義用到Days,所以定義在前面
set D := {1..Days};
set M := {20, 50, 100, 180, 300, 540, 780, 1020};
set S := {9..12};
set C := {"point", "paper", "candy", "bag", "moji"};
param Max_ap := 1800;
對於數學表記當中帶有「下標」的參數或變數,ZIMPL的宣告方式是在名稱後面接上[],裡頭放著一個集合。比如對每個$s \in S$我們有$Ap\_cost_s$,在ZIMPL宣告時會寫成Ap_cost[S]。那要怎麼指定個別的值呢?直接寫出來可能比較好懂:
param Ap_cost[S] := <9> 20 <10> 20 <11> 20 <12> 20;
如果沒有指定每一個值的話,則可以寫default xx代表「沒寫出來的值都是xx」。像上面這個例子可以直接寫default 20。
對於下標有兩種代號的,比如$Stage\_get_{s, c}$,我們得將$s, c$視為一個從集合$S \times C$來的元素(笛卡兒積),而在ZIMPL裡頭這個集合會用乘號寫成S * C。
要指定值的時候,除了用上面的表達方法,將一個下標寫成<s, c>之外,也可以用二維表格來呈現:
param Stage_get[S * C] :=
| "point", "paper", "candy", "bag" |
| 9| 32 , 65 , 0 , 0 |
| 10| 32 , 0 , 64 , 0 |
| 11| 32 , 0 , 0 , 46 |
| 12| 32 , 23 , 25 , 19 | default 0;
可以看到$S$的部分寫在左側,$C$的部分寫在上面,可以比較直觀看出諸如$Stage\_get_{9, point} = 32$的值。
在這個例子沒有直接指定的值,比如$Stage\_get_{9, moji}$,則會被設定成default的0。
這邊已經把參數設定時需要知道的文法講過一次了,接下來就把剩下的參數也指定完畢:
param Need[C] :=
<"point"> 35005,
<"paper"> 17350,
<"candy"> 14277,
<"bag"> 11568,
<"moji"> 226;
param Rate_paper2candy := 5;
param Rate_candy2bag := 5;
param Rate_bag2moji := 200;
param M_get[C] :=
<"point"> 120,
<"paper"> 40,
<"candy"> 35,
<"bag"> 30,
<"moji"> 1;
由於Python是主流的程式語言,筆者不會對語法有太多著墨,更多的是針對PySCIPOpt的用法進行介紹。
首先要先創造一個Model物件:
from pyscipopt import Model
model = Model("我要拿5星フブキ的線性規劃問題") #名字隨便取
接著宣告集合與參數。
這些參數會是普通的Python變數(不要搞混Python變數與先前提到的「線性規劃的變數」),因為對解題器來說這些值是固定的,我們在之後寫限制式時直接調用Python變數的值即可。
雖然說是「集合」但我們用了list去儲存,算是一種便宜行事,因為對Python使用者來說list還是比較熟悉的存在。我們使用的時候也只需要知道集合內「有哪些」元素即可,只要這裡頭的元素真的沒有重複,那並不會發生問題。
Days = 13
D = list(range(1, Days+1))
M = [20, 50, 100, 180, 300, 540, 780, 1020]
S = list(range(9, 12+1))
C = ["point", "paper", "candy", "bag", "moji"]
Max_ap = 1800
對於數學表記當中帶有「下標」的參數或變數,我們可以用dict儲存,並用comprehension的技巧一口氣定義所有的鍵值與預設值:
Ap_cost = {s: 20 for s in S}
Python也內建笛卡兒積product:
from itertools import product
Stage_get = {sc: 0 for sc in product(S, C)}
Stage_get[9, "point"] = 32
Stage_get[9, "paper"] = 65
Stage_get[10, "point"] = 32
Stage_get[10, "candy"] = 64
Stage_get[11, "point"] = 32
Stage_get[11, "bag"] = 46
Stage_get[12, "point"] = 32
Stage_get[12, "paper"] = 23
Stage_get[12, "candy"] = 25
Stage_get[12, "bag"] = 19
有時我們也不用comprehension,比如在每個值都不一樣的狀況:
Need = {
"point": 35005,
"paper": 17350,
"candy": 14277,
"bag": 11568,
"moji": 226
}
剩下的參數:
Rate_paper2candy = 5
Rate_candy2bag = 5
Rate_bag2moji = 200
M_get = {
"point": 120,
"paper": 40,
"candy": 35,
"bag": 30,
"moji": 1
}
輸入變數
宣告變數的關鍵字是var。
對於整數型態的變數我們需要事先表明它是整數integer,我們也可以在宣告時直接指定變數的基本範圍:比如說非負整數就是integer >= 0。
var ap[D] integer >= 0;
var time[D * S] integer >= 0;
對於值域為0或1的整數變數,有個專有名詞叫做binary:
var to_ap[D * M] binary;
有關直接指定binary、指定integer >=0 <=1或者用限制式限制的差異
在宣告變數時寫出integer是告訴解題器這個變數的性質是整數,這代表它不需要考慮實數的稠密性之類的性質,只需要考慮整數的性質即可。定義範圍則是告訴解題器可以直接只考慮某個特定範圍,與限制式的差異是:在解題器還在處理其它限制式時,對這變數的假設就會直接以這個比較小的範圍做考量,而不是在看到限制式時才限縮範圍。
那integer >=0 <=1與binary又有什麼差別?在做整數規劃(Mixed Integer Programming)時我們常常使用到各種0或1的變數,這些變數同時擔任「開關」的任務,在某些表述語言或解題器甚至會利用這種開關變數來實作「條件分歧」的語法糖。
如果寫成integer >=0 <=1的形式,不是每個地方都會貼心地轉成開關變數,更別說如果寫成「與1等價的某個複雜式子」時還要看出來這個式子與1等價……總之要用的時候寫binary就對了。
剩下的變數:
var base_get[C] integer >= 0;
var ext_get[C] integer >= 0;
var total_get[C] integer >= 0;
var dummy[C] integer >= 0;
接下來我們要來宣告餵給SCIP的「變數」了。
對模型加變數的方法就是呼叫addVar,它同時也會回傳一個Python變數,除了讓我們可以拿著這個Python變數去跟模型索取對應的「模型變數」的資訊之外,更重要的是在之後列式時我們也可以直接拿著這些Python變數做Python數學運算,模擬模型內的數學運算。
稍微解釋一下addVar的用法,一般常用到addVar的前四個參數:addVar(name = '', vtype = 'C', lb = 0.0, ub = None)name代表名稱,這個可以隨便取,會決定在模型輸出答案時顯示的變數名稱vtype代表型態,'C'為連續(實數)、'I'為整數,'B'則為0或1的整數lb為下限,如果設定為None則為負無窮大ub為上限,如果設定為None則為正無窮大
ap = {model.addVar(f"ap_{d}", 'I') for d in D}
因為預設的下限就是0、無上限,所以後面兩個參數都不用設定。
其它變數也能設定:
time = {ds: model.addVar(f"time_{ds}", 'I') for ds in product(D, S)}
to_ap = {dm: model.addVar(f"to_ap_{dm}", 'B') for dm in product(D, M)}
base_get = {c: model.addVar(f"base_get_{c}", 'I') for c in C}
ext_get = {c: model.addVar(f"ext_get_{c}", 'I') for c in C}
total_get = {c: model.addVar(f"total_get_{c}", 'I') for c in C}
dummy = {c: model.addVar(f"dummy_{c}", 'I') for c in C}
對於0或1的變數,指定型態為'B'就可以使之用在一些地方簡化式子(後述)。
目標式
宣告目標式的關鍵字是minimize或maximize。
儘管我們沒有為我們的式子取代號,但在ZIMPL我們需要幫式子取名字。
取完名稱後要加上冒號。
接著我們看看如何表達$\sum_{d \in D} ap_d$這個算式。
$\Sigma$寫成sum,後面接它要看的集合與元素名稱:<d> in D。
之後接冒號,後面接$ap_d$的部分。
先前宣告$ap_d (d \in D)$的時候我們是寫var ap[D],在中括號中間夾集合名稱。而要表達特定一個$ap_d$,就是把中括號裡頭的東西改成集合內的個別元素d。
minimize cost: sum <d> in D: ap[d];
對模型加上目標式的方法是呼叫setObjective。它只做設定並不回傳值。
一般常用到的前四個參數:setObjectivesetObjective(coeffs, sense = 'minimize')代表一個「係數(coefficient)」coeffs代表這個目標是senseminimize或是maximize
這邊所謂的「係數」可以是一個變數,也可以是包含加減乘除等運算的算式。也就是運算結果是一個「數值」的式(而非「真假值」)。
(如果先前的數學列式都是正確的,那麼這個小細節不會害你出事。)
接著我們看看如何表達$\sum_{d \in D} ap_d$。
$d \in D$的部分,用comprehension表達即可。
Python雖然已經內建sum函數可以讓我們方便寫出$ap_1+ap_2+…$,但PySCIPOpt建議我們使用它定義的quicksum代替sum的使用,原因是這樣子SCIP看到的就是一個$\Sigma$,而非拆開來的$ap_1+ap_2+…$,代表SCIP為$\Sigma$進行的優化得以進行。
from pyscipopt import quicksum
model.setObjective(quicksum(ap[d] for d in D), 'minimize')
限制式
宣告限制式的關鍵字是subto。同樣地我們也需要幫限制式取名。語法跟目標式基本相同。
要表達$\forall d \in D.$,需要寫成forall <d> in D:。
$=$寫成==、$\le$寫成<=,$\ge$則是>=。
先前提到如何表達$ap_d$,現在要表達$time_{d, s}$。當時用time[D * S]定義,現在中括號裡頭的兩個元素d與s要用逗號分隔,寫成time[d, s]。
subto ap_cost:
forall <d> in D:
ap[d] == (sum <s> in S: Ap_cost[s] * time[d, s]);
subto everyday:
forall <d> in D:
ap[d] <= Max_ap;
針對$\forall d \in D, m \in M. ap_d \ge m \times to\_ap_{d, m}$這條式子,我們當然可以直接寫:
subto ap_atleast:
forall <d, m> in D * M:
ap[d] >= m * to_ap[d, m];
不過在我們列這條數學式時曾經提過,「如果這個開關為1則某件事情要成立」是做線性規劃時很常見的一種情境。在這邊,「如果$to\_ap_{d, m} = 1$則$ap_d \ge m$」與「如果$to\_ap_{d, m} = 0$則$ap_d \ge 0$」。後者實際上就是「沒有額外限制」,因為$ap_d$本來就是大於等於0。
這種「只有當變數處於某種情形出現時才需要成立的條件式」,在ZIMPL有個方便的語法vif-then[-else]-end,比如上式即可改寫成:
subto ap_atleast:
forall <d, m> in D * M:
vif to_ap[d, m] == 1 then ap[d] >= m end;
有關vif使用的注意事項
vif後面接的條件式只能包含整數變數(integer與binary)的線性組合(可以包含沒有等號的$< \ne >$以及所有種類的邏輯運算子),並且出現的整數變數必須在宣告時就限制在有限的範圍內。意思是參與的變數不能是:
var foo integer >= 0;
subto foo_range: foo <= 5;
而必須是:
var foo integer >= 0 <= 5;
這個變數才能用在vif。binary變數本身就只有兩種可能性,所以沒這個問題。
在之前講到變數$dummy$時曾經帶過這個話題,現在這邊可以一窺兩者的差異。不過我們並沒有要在vif後面接$dummy$,所以我們不會碰上這個議題。
至此我們已經把我們會用到的語法解釋完了。剩下的限制式:
subto total:
forall <c> in C:
base_get[c] + ext_get[c] == total_get[c];
subto need:
forall <c> in C:
total_get[c] >= Need[c];
subto base:
forall <c> in C:
base_get[c] ==
(sum <s> in S:
(sum <d> in D: time[d, s]) * Stage_get[s, c]
) +
(sum <d, m> in D * M: to_ap[d, m]) * M_get[c];
subto ext_point: ext_get["point"] == 0;
subto ext_paper: ext_get["paper"] == 0;
subto ext_candy: ext_get["candy"] == (total_get["paper"] - Need["paper"] - dummy["paper"]) / Rate_paper2candy;
subto dummy_paper: dummy["paper"] <= Rate_paper2candy - 1;
subto ext_bag: ext_get["bag"] == (total_get["candy"] - Need["candy"] - dummy["candy"]) / Rate_candy2bag;
subto dummy_candy: dummy["candy"] <= Rate_candy2bag - 1;
subto ext_moji: ext_get["moji"] == (total_get["bag"] - Need["bag"] - dummy["bag"]) / Rate_bag2moji;
subto dummy_bag: dummy["bag"] <= Rate_bag2moji - 1;
添加限制式的方法是呼叫addCons。它雖然同時也會回傳一個指向模型內部的限制式的參照,不過我們一般用不太到,姑且可以無視其回傳值。
addCons的用法,一般只需要知道前兩個參數:addCons(cons, name = '')
注意,與addVar不同,名稱是擺在第二個位置cons代表一個「限制(constraint)」
這邊所謂的「限制」是由>=、==與<=與左右兩邊各一個「係數」組合而成。
(如果先前的數學列式都是正確的,那麼我們的式子不會出問題。)
在數學式我們會用$\forall$代表許多限制式,而在PySCIPOpt用for迴圈逐一添加是最直覺的:
for d in D:
model.addCons(ap[d] == quicksum(Ap_cost[s] * time[d, s] for s in S), f"ap_cost_{d}")
for d in D:
model.addCons(ap[d] <= Max_ap, f"everyday_{d}")
針對$\forall d \in D, m \in M. ap_d \ge m \times to\_ap_{d, m}$這條式子,我們當然可以直接寫:
for d, m in product(D, M):
model.addCons(ap[d] >= m * to_ap[d, m], f"ap_atleast_{d}_{m}")
在我們列這條數學式時曾經提過,「如果這個開關為1則某件事情要成立」是做線性規劃時很常見的一種情境。在這邊,「如果$to\_ap_{d, m} = 1$則$ap_d \ge m$」與「如果$to\_ap_{d, m} = 0$則$ap_d \ge 0$」。後者實際上就是「沒有額外限制」,因為$ap_d$本來就是大於等於0。
跟ZIMPL比較自由的前提寫法不同,PySCIPOpt的寫法限定前提是「開關為1(或0)時」的形式:addConsIndicator(cons, binvar = None, activeone = True, name = "IndicatorCons")cons是前提成立後要成立的條件,不能是$a \le foo \le b$這種同時具有上限與下限的條件(若真的有需求則必須分開寫)binvar是0或1的變數,宣告時需要標示vtype = 'B'。如果None則內部會自動創造一個二元變數activeone如果為True則前提是「若開關為1時」,False則為「若開關為0時」
以基本的線性規劃而言基本上掌握「若開關為1/0時則……」應該就夠用了。
for d, m in product(D, M):
model.addConsIndicator(ap[d] >= m, to_ap[d, m], name = f"ap_atleast_{d}_{m}")
剩下的限制式基本上就沒有什麼要注意的地方了:
for c in C:
model.addCons(base_get[c] + ext_get[c] == total_get[c], f"total_{c}")
for c in C:
model.addCons(total_get[c] >= Need[c], f"need_{c}")
for c in C:
model.addCons(base_get[c] ==
quicksum(quicksum(time[d, s] for d in D) * Stage_get[s, c] for s in S) +
quicksum(to_ap[d, m] for d, m in product(D, M)) * M_get[c],
f"base_{c}")
model.addCons(ext_get["point"] == 0, "ext_point")
model.addCons(ext_get["paper"] == 0, "ext_paper")
model.addCons(ext_get["candy"] == (total_get["paper"] - Need["paper"] - dummy["paper"]) / Rate_paper2candy, "ext_candy")
model.addCons(dummy["paper"] <= Rate_paper2candy - 1, "dummy_paper")
model.addCons(ext_get["bag"] == (total_get["candy"] - Need["candy"] - dummy["candy"]) / Rate_candy2bag, "ext_bag")
model.addCons(dummy["candy"] <= Rate_candy2bag - 1, "dummy_candy")
model.addCons(ext_get["moji"] == (total_get["bag"] - Need["bag"] - dummy["bag"]) / Rate_bag2moji, "ext_moji")
model.addCons(dummy["bag"] <= Rate_bag2moji - 1, "dummy_bag")
答案解析
將撰寫好的模型儲存成.zpl檔,並且呼叫SCIP:
scip -f <檔案名稱>.zpl
將會輸出一大堆資訊。我們來找最重要的部分。我們首先找到SCIP Status的部分:
SCIP Status : problem is solved [optimal solution found]
Solving Time (sec) : 0.00
Solving Nodes : 4
Primal Bound : +2.03400000000000e+04 (2 solutions)
Dual Bound : +2.03400000000000e+04
Gap : 0.00 %
只要顯示的是problem is solved [optimal solution found]就代表有找到一組解答。
如果是problem is solved [infeasible]則代表找不到符合條件的解,如果顯示這個,在這題通常是代表$Max\_ap$設定太低了,代表沒有辦法在每天只花那麼少AP的情況下達成所有條件。
回到這組解答。我們再關注Primal Bound的部分,顯示的是20340($2.0340 \times 10^4$),代表我們的目標式$\sum_{d \in D} ap_d = 20340$。
接著要分析我們之前定義的各種變數。以我們實際刷關卡的角度來看,我們最關注的莫非:要打關卡幾多少次,也就是$time$變數的部分。
從SCIP Status的部分向下找到primal solution (original space)的部分會列出所有由我們定義的變數的值。雖然顯示得有些雜亂,但抓出所有名稱為time開頭的變數應該不算太難:
time#10#12 51 (obj:0)
time#11#12 51 (obj:0)
time#12#12 58 (obj:0)
time#13#9 51 (obj:0)
time#13#12 39 (obj:0)
time#1#11 90 (obj:0)
time#2#11 90 (obj:0)
time#3#11 90 (obj:0)
time#4#11 90 (obj:0)
time#5#11 90 (obj:0)
time#6#12 90 (obj:0)
time#7#12 90 (obj:0)
time#8#11 86 (obj:0)
time#9#12 51 (obj:0)
這組解告訴我們「第一天刷90次關卡11」、「第10天刷51次關卡12」等資訊。
但我們知道,51或90的重點其實只在於每天都至少刷滿51次可以觸發所有突襲事件,以及每天花費的AP數量被限制在90場1800,詳細第幾天到底刷51場還是90場其實並沒有那麼重要,把第一天與第10天的結果互換只不過是變成「另一組解」罷了。
這些資訊更寶貴的地方是告訴我們:「在這13天內,我們總共需要刷51次關卡9、536次關卡11與430次關卡12」(這組解說不用打關卡10),刷這麼多次就可以讓我們拿到需要的所有通貨量。我們只需要在每天刷1020AP的前提下在13天內刷滿這目標即可。
小技巧
要求出打某個關卡的總次數,可以加上額外的變數$total\_time_s = \sum_{d \in D} time_{d, s}$來監控:
var total_time[S] integer >= 0;
subto total_times:
forall <s> in S:
total_time[s] == sum <d> in D: time[d, s];
這樣一來看$total\_time_s$即可馬上得知關卡$s$的總次數。
但是這個技巧需要注意觸發突襲事件的次數。筆者所使用的參數使得我必須每天刷滿1020AP,但如果需求量沒有那麼大,可能有些日子不需要刷滿8場突襲,則需要自行調整刷關計畫。
模型執行此方法:
model.optimize()
將會輸出一部分資訊並將詳細的資訊儲存於內部。
首先我們檢查有沒有解,觀察剛剛的輸出,下方會有一行:
SCIP Status : problem is solved [optimal solution found]
或者
model.getStatus()
如果是optimal則代表找到一組解。
如果是problem is solved [infeasible]或getStatus()結果為infeasible則代表無解。如果是這個狀況,在這題通常是代表$Max\_ap$設定太低了,代表沒有辦法在每天只花那麼少AP的情況下達成所有條件。
回到這組解答。我們再關注Primal Bound的部分,顯示的是20340($2.0340 \times 10^4$),代表我們的目標式$\sum_{d \in D} ap_d = 20340$。
我們也可以用:
model.getObjVal()
得到這個數字。
接下來是分析這組解裡頭各個變數的資訊。首先我們把現在這組解先存下來:
sol = model.getBestSol()
要查看各個變數的值,只需要呼叫:
model.getSolVal(sol, ap[1])
這樣就可以查看ap[1]的值。
要查看每天花費多少AP:
for d in D:
print(f"ap_{d} = {model.getSolVal(sol, ap[d])}")
輸出:
ap_1 = 1800.0
ap_2 = 1800.0
ap_3 = 1020.0
ap_4 = 1800.0
ap_5 = 1020.0
ap_6 = 1800.0
ap_7 = 1040.0
ap_8 = 1800.0
ap_9 = 1800.0
ap_10 = 1800.0
ap_11 = 1800.0
ap_12 = 1800.0
ap_13 = 1060.0
可以看出每一天都要至少花1020AP,代表每天8場突襲都是打好打滿的。
以我們實際刷關卡的角度來看,我們最關注的莫非:要打關卡幾多少次,也就是$time$變數的部分。所以可以運行:
for d in D:
for s in S:
print(f"time_{d},{s} = {model.getSolVal(sol, time[d, s])}")
去查看各個變數的值。
但我們已經知道每天8場突襲都要打滿,每天第幾關各打幾場的重點其實只在於每天都至少刷滿51次可以觸發所有突襲事件,以及每天花費的AP數量被限制在90場1800,詳細第幾天到底刷51場還是90場其實並沒有那麼重要,把某兩天的結果互換只不過是變成「另一組解」罷了。
我們再試著印出關卡次數的加總:
for s in S:
print(f"total_time_{s} = {sum(model.getSolVal(sol, time[d, s]) for d in D)}")
輸出是:
total_time_9 = 53.0
total_time_10 = 1.0
total_time_11 = 539.0
total_time_12 = 424.0
也就是說這組資訊最寶貴的地方是告訴我們:「在這13天內,我們總共需要刷53次關卡9、1次關卡10、539次關卡11與424次關卡12」,刷這麼多次就可以讓我們拿到需要的所有通貨量。我們只需要在每天刷至少1020AP的前提下在13天內刷滿這目標即可。
如果你不是刷關卡9~12,而是8~11之類,AP消耗數有所不同的關卡
上文假設我們刷的是關卡9~12,每個關卡都是消耗20AP。但如果新手打不過關卡12只能退守8~11,或關卡11也不行,只能刷7~10呢?
這樣在特定狀況下會有個細節需要注意,怎麼說呢?
假如只有2天,然後規劃出來是每天各刷1020AP,但第一天是關卡8(15AP)刷68次,第二天是關卡9刷51次。我們說要看刷關卡的總次數與突襲總次數,而這裡我們不能「第一天刷67次關卡8與1次關卡9(1025AP)、第二天刷1次關卡8與50次關卡9(1015AP)」,這樣儘管刷關總次數是對的,但突襲總次數錯了。
一般來說一組題目的解答不會只有一種(又或者找不到解答)。但每組解都確實滿足我們列出的所有條件,所以我們只需要關注一組解即可。如果看完解答之後很在乎某件事情,那麼可以自行更改參數或增加條件式並且重新計算解答,得到自己最滿意的結果。
後記
本網誌先前設置於NexusBytes的虛擬伺服器中,同標題的「本文」也於1月春節假期之前完稿並發表,但於發表後不久伺服器即無預警關閉,求助客服無門,應該是人間蒸發了(當初官網仍然持續運作,繼續販賣著已經不再運作的服務,目前網站已經消失)。又因備份間隔過長導致文章沒有成功備份,導致我搬移伺服器之後文章必須重寫。在搬移到位於ABLENET的新伺服器之後,重新思考文章架構,才有了現在版本的這篇文章。如果有疏漏之處亦會隨時編輯補充。