一種解決藥店賣藥問題的新思路

期末在做數位鍾，在 csdn 上參考了前人的經驗。遂想到把資料結構大作業的實驗報告也發上來，興許可以供學弟學妹參考。

解決這道藥店賣藥問題最普遍也是最容易想到的做法是貪心求解，然而本題多限制之間相互作用，關係複雜，因此簡單的貪心難以找到十分優秀的解。本文提供了一種不同的思路，從新的角度對這個問題進行思考，設計了新的演演算法，並且詳述了演演算法主要思路產生過程及其合理性。實際執行結果優秀，亂資料下（也許大概率）能找到全域性最優解。

報告注重思路的分析和證明，實現部分稍簡略，具體可見程式碼。

1 實驗背景及要求

你是⼀家藥店的老闆，這個⽉你從供貨商⼿⾥收到了⼀批共50個藥品，其中每個藥品有獨立的進價和保質期，其中進價的區間為[20,30]元，保質期的剩餘天數為[1,15]天.你每天可以陳列出10個藥品在商店中售賣,每天會有三個顧客各過來購買⼀個藥品。藥品的定價只能為商品進價加上{-1.5,-1, -0.5, 0, 2 ,4 ,6}元，不得隨意定價。每位顧客購買的邏輯是，買其中最便宜的藥品，如果說最便宜的藥品價格⼀致則購買保質期⻓的藥品。三位顧客會依次購買藥品。

藥品如果沒有陳列在商店中，⽽是存放在倉庫時，會收取管理費，其中保質期低於5天的每天收取1元管理費，其餘的每天收取0.5元。每天的陳列藥品可以少於10個

你的⽬標是，10天之後，你的利潤（售出商品售價總和-售出商品進價總和-⽀付倉庫的管理費⽤-10天內過期/丟棄商品的進價）最⼤。

2 任務分析

2.1 初步認識

首先，肯定不能通過直接列舉的方法求解，因為如果把每一種藥品擺放方式及定價策略看做一個解，那麼藥品選擇方案數為 ${\prod }_{i=0}^9\left(\genfrac{}{}{0ex}{}{50-3i}{10}\right)$，定價方案數為 $6^{100}$，因此解空間大小為 ${\prod }_{i=0}^9\left(\genfrac{}{}{0ex}{}{50-3i}{10}\right)\cdot 6^{100}$，巨大無比。一種解決思路是貪心，想出一些看上去靠譜的策略，然後根據這些策略做唯一選擇。這樣做的好處是效率較高，程式碼容易實現；壞處是難以考慮所有情況，不能保證結果足夠優秀。

我的思路是：首先觀察問題性質，根據其性質縮小解空間，後對解空間進行搜尋。

不管怎麼說，要解決這個問題，就要對其有一個宏觀的認識。注意到存在以下幾個平凡的結論：

1. 每天都會賣出三種藥。
2. 我們傾向於先賣即將到期的那些藥，因為他們會產生額外的管理費。
3. 有時候會刻意擺不滿貨架，來引導顧客選擇我們想賣給他們的藥。
4. 我們傾向於先賣便宜的藥，因為如果要先賣貴的藥，要麼就通過降低售價使之被選走，要麼就要減少擺上貨架的藥品數量增加管理費。而先賣便宜的就沒有這麼多顧慮。
5. 如果一種藥最後會過期或者被扔掉，那麼一定是在第一天的開始就扔掉。

基於以上觀察，能想到幾種貪心策略，比如最簡單的把所有藥品按照過期時間為第一關鍵字、價格為第二關鍵字升序排序，直接按照這個順序選擇藥品，並把定價都設為最高。這樣確實能夠得到一組不錯的解，但我們的追求遠不止於此。

2.2 進一步思考

經過分析，我們可以接著得出以下性質：

性質2.2.1：最優方案一定是扔掉的藥數量最少的某一種方案。

證明：扔掉一種藥至少虧損20元成本價。記 $f_n$ 為扔 $n$ 種藥的最大獲利，我們採用反證法，假設有 $f_k<f_{k+1}$，那麼 $f_{k+1}$ 對應的賣藥方案中肯定存在一天可以賣最後被扔掉的藥但是沒賣。現在我們在這一天的貨架上只擺出想要賣的三種藥，其他位置為空，這樣當天多賺了至少二十塊成本，虧了至多7*1元的管理費，同時因為不扔藥最多額外產生7.5元管理費。因此總盈利 $profit\ge f_{k+1}+20-7\times 1+5\times 1+5\times 0.5>f_{k+1}$，與 $f_{k+1}$ 為最大獲利矛盾。因此，$f_k\ge f_{k+1}$ 恆成立，證畢。

性質2.2.2：如果一種藥 $(v_i,t_i)$，另一種藥 $(v_j,t_j)$，滿足 $v_i<v_j$ 並且 $t_i<t_j$，那麼把 $i$ 放在 $j$ 前面肯定不會比把 $j$ 放在 $i$ 前面更劣。

證明：如果過期時間相同，肯定是先賣便宜的；如果價格相同，肯定先賣過期早的。

根據性質2.2.1，我們要儘量把將要過期的藥都賣掉。首先把藥按照過期時間升序排序，記 $s_{l,r}$ 表示過期時間在第 $l$ 天到第 $r$ 天之間的藥品集合。注意到，如果 $|s_{1,x}|<3x$，最優策略肯定是將 $s_{1,x}$ 中的藥品全部賣出；否則，過期時間在 $[1,x]$ 的藥肯定不能都售出，並且從第 1 天到第 $x$ 天賣的一定都是 $s_{1,x}$ 中的藥，此時第 1 天到第 $x$ 天賣的藥和別的時間賣的藥是相對獨立的，分開計算並相加即可。

進一步地，如果 $|s_{1,x}|\ge 3x$ 並且已經處理完了 $[1,x]$ 天，那麼接著把第 $x+1$ 天作為第一天重複上述過程即可。

推論2.2.1：可以把排序後的藥品劃分為若干段，每段的選擇相互獨立。

過期時間隨機生成的性質使得這個「分裂」的過程大概率會進行，從而縮小問題規模。

2.3 壓縮解空間——優雅的列舉

再次從一個宏觀的視角思考，該如何確定我們的策略？總的來說，為了縮小解空間，我們要對關鍵資訊，即對結果影響大的變數進行決策。做一個不恰當的類比，就是要找到線性空間的一組基。一般來說有以下兩種思路：

1. 先確定每一種藥品的定價，再確定每天擺什麼藥。
2. 先確定每天擺什麼藥，再給要定價。

而對於每天擺什麼藥，還有兩種思路：

1. 直接確定擺哪些藥，根據定價即可得知賣出了哪三種藥。
2. 先確定賣出哪三種藥，再決定擺哪些藥。

容易發現，先確定每一種藥品定價是不怎麼優秀的，因為最後總共只會賣出 30 瓶藥，別的藥的價格無關緊要（擺上貨架不賣的藥肯定越貴越好）；而且一種藥在不同時間定價可以不同，難以預先確定。

而對於擺藥，沒賣出去的擺上貨架的藥其實沒有那麼重要。假設每天賣哪三種藥及其定價已經確定了，那麼只要再擺一些售價高於它們的藥即可（當然優先擺五天內過期的）。而對於定價，只需要列舉賣出的三種藥的定價，即可算出當天收益，取最優即可。

總結一下整體思路：首先確定每天賣哪三種藥，然後逐天給這些藥定價，順便確定擺那些藥能剩下最多的管理費。這樣做的優勢是，我們只需要列舉每天賣哪三種藥和這三種藥的售價，其他的量是可以直接計算出來的。實際上我們只需要列舉賣出的三種藥裡的最高售價即可。

再進一步，根據性質1.2.2，只需要列舉賣的是哪一天過期的藥。粗略估算解空間一個很鬆的上界為 ${\left(\genfrac{}{}{0ex}{}{15}{3}\right)}^{10}$，遠遠小於原本的大小，而且其中仍有大量的違反性質2.2.1不優的解，可以被直接丟棄。最為關鍵的是這種問題規模的壓縮是「無失真」的，因為扔掉的都是那些不優的解，而最優解，或者說比較優的那些解一定還在剩下的解空間中。

事實上，如果過期時間在 [1, 15] 內均勻分佈，進一步用性質2.2.1和推論2.2.1削減後的解空間大小小於 1e8 的概率約為 70%（實測），也就是說直接列舉就可以輕鬆找到全域性最優解。可惜的是，過期時間並不是完全隨機的，因此需要另尋他法。

2.4 搜尋解空間——模擬退火

如上面分析，對於解空間的縮小已經做到我能力範圍內的極限了，仍然不可能通過窮舉解決，因此不得不採用隨機演演算法對其進行搜尋。我採用了模擬退火的方法。

具體來說，將一張3行10列且每列無序的數表看做一個解，先按照每天賣出最早過期的三種藥品的方式貪心得到初始解。每次調整時，有兩種方法：1. 隨機某一天的某一種藥，和另一天的一種藥交換；2. 隨機某一天的某一種藥，和某一種沒被選擇的藥交換。注意，一定滿足調整後的序列滿足性質2.2.1，即不能產生新的過期的藥，否則會有大量無效搜尋。做如上調整 1~5 次，生成新的解。

3 任務實現

3.1 模擬部分

3.1.1 關鍵變數

函數/變數名稱	型別	功能
read()	function	解析命令列引數並讀取相應檔案
simulate()	function	模擬賣藥過程
a	node*	藥品結構體，存放價格(v)和過期時間(t)
sold	int*	是否被賣出/扔掉
del	vector	每天扔掉哪些藥

3.1.2 實現流程

寫在 sim.cpp 裡。逐天處理，先扣除每個藥的管理費，再挨個上架藥品，上架的同時把扣掉的管理費加回來。排序得到賣出的三件藥品，計算收益後將 sold 設為 1。最後把當天應該扔掉的藥品扔掉，扣除相應進價。

3.2 制定策略

3.2.1 關鍵變數

函數/變數名稱	型別	功能
SA()	function	模擬退火流程
calculate()	function	計算某種策略的盈利
init()	function	貪心生成模擬退火初始解
adjust()	function	從某個解隨機調整生成新的解
R	Random	亂數類
state	SAnode	模擬退火中用到的策略類（每天賣哪三種）
ksj1, ksj2	vector	計算盈利時分別存放過期時間<=5和>5的藥

3.2.2 實現流程

簡述模擬退火流程：

1. 生成初始解 state，設定初始溫度 T
2. 隨機調整解，計算新解的盈利
3. 如果新解更優，就接受這個解；否則按照 $P=e^{\frac{\Delta profit}{T}}$ 的概率接受這個解
4. 降溫，若達到最低溫度退出，否則回到步驟2

關於退火引數的設定本人沒有什麼經驗，完全憑感覺亂設的。但我猜測最優解距離初始貪心解不會太遠，因此初始溫度設定不高。此外，在最後我還對所得最優解的附近進行爬山，企圖找到一組更優的解。

這個程式不需要用到上面的 sim.cpp，而是在 calculate() 中確定每天藥品擺放和定價的時候順便求出收益，因此此部分寫在 s1.cpp 中。

4 實驗資料

資料	data1	data2	data3	data4	data5
教輔	24	27.5	22.5	27.5	23.5
本程式	36.5	35	30.5	39.5	31

5 缺點和不足

在扔掉那些必須要扔的藥的時候，我是選擇儘量便宜的幾種藥丟棄，這樣做在絕大多數情況下都是很好的。但在構造資料下，可能會出現扔掉較貴的藥反而獲利更大的情況。這也是我的程式中的唯一紕漏，也是其將最優解排除在解空間外的唯一情況。綜合考慮亂資料下這種情況微乎其微的概率以及完善它所需要的程式碼量，沒有選擇去修補它。

6 總結與感悟

從結果看，執行多次程式結果相同，這令我非常吃驚，因為我的程式裡有隨機化成分，在我預想中最終結果應該有微小浮動，至少不會總是一樣。對此有兩種可能的解釋：一種解釋是這個解就是全域性最優解，自然不會有更好的解出現；另一種解釋是這是區域性最優解而全域性最優解極為刁鑽，甚至不能近其身。我個人傾向於第一種。因為每多一種藥過期在我的解法裡相當於多了一條限制，會幫助我縮小解空間。而一開始下發的data1中沒有藥過期，理論上這種情況下我的演演算法找到最優解的概率最小，但它僅用很短的時間就找到了一個解，並且通過手算可以證明它就是最優解。

當然，這也無法完全驗證或者證偽，因為目前還沒有別人跑出更優的解。究其原因，可能是2.3中給出的解空間上界確實非常非常鬆，實際資料下解空間大小遠沒有這麼大，從而比較容易找到最優解；也有可能是最優解距離初始貪心解的距離不太遠。

回過頭來看，整個演演算法最為關鍵的部分是前面對於壓縮解空間的分析，後面的模擬退火反而沒那麼重要，只是順勢而為罷了，或者只能算一種妥協。事實上，我把模擬退火的部分去掉後，直接從初始貪心解周圍隨機若干解的效果也不錯，甚至也超過了教輔的收益。將隨機換成直接爬山效果更佳，距離模擬退火退出來的結果只有幾元錢的差距。

最終的 sim.cpp 和 s1.cpp 加上實驗報告總共花了三個晚上的時間，比預期短，是因為我吸取了銀行的教訓，沒有過於自信直接上手而花費很多時間在不可行的方案上。這次我在前期花了大量的時間準備思考，寫了大量測試資料的程式碼（雖然這部分最後沒派上用場），全部考慮清楚後再著手寫程式碼，避免了很多無用功，也讓我的思路更加清楚，程式碼一氣呵成，沒花多少時間偵錯。

/**
 * TODO:
 * 
 * 
 * */

# include <bits/stdc++.h>
# include <getopt.h>
# define pb push_back
# define ll long long
using namespace std;
const double deltav[] = {-1.5, -1, -0.5, 0, 2, 4, 6};  // 定價表
struct node {
    double v;
    int t;
} a[100];

string fpath = "testdata\\data";  // 路徑，測試用

// 策略結構體，包含strategy和delete
struct Koutput {
    int s[11][11][2], d[51][2], n;
    Koutput() { n = 0; }
    Koutput& operator=(const Koutput &b) {
        memcpy(s, b.s, sizeof(s));
        memcpy(d, b.d, sizeof(d));
        n = b.n;
        return *this;
    }

    void print() {
        ofstream fout;
        fout.open((fpath + "\\strategy.my").c_str());
        for (int i = 1; i <= 10; i++) {
            for (int j = 0; j < 10; j++) {
                fout << s[i][j][0] << ',' << s[i][j][1] << ' ';
            }
            fout << '\n';
        }
        fout.close();
        fout.open((fpath + "\\delete.my").c_str());
        for (int i = 0; i < n; i++) {
            fout << d[i][0] << ' ' << d[i][1] <<'\n';
        }
        fout.close();
    }
} beststr;

// 退火用到的3*10的數表結構體
struct SAnode{
    int A[11][3], B[30], n;
    SAnode(): n(0) {}
} initstate, bestSAstr;

int n = 50;  // 藥品總數
int s_a[100];  // 過期時間升的藥品編號
int initDel[100], initDelcnt;  // 一開始扔掉哪些藥及其數量
double bestans = -10000;  // 最大獲利
double initCost = 0;  // 扔藥的損失
vector<int> ksj1, ksj2;  // 計算獲利時維護的集合，ksj1是即將過期的藥，ksj2是其餘藥

bool cmp_greater(int x, int y) { return a[x].v > a[y].v; }
bool cmp_id(int x, int y) {
    if (a[x].t ^ a[y].t) return a[x].t < a[y].t;
    return a[x].v > a[y].v;
}

// 生成[0,1]浮點數
double frand() { return 1.0 * rand() / RAND_MAX; }

// 解析命令列引數並讀入
void read(int argc, char* argv[]) {
    if (argc == 1) {
        int tmp;
        scanf("%d", &tmp);
        fpath = fpath + to_string(tmp);
        FILE* fp = fopen((fpath + "\\medicine.txt").c_str(), "r");
        for (int i = 0; i < n; i++) {
            fscanf(fp, "%lf%d", &a[i].v, &a[i].t);
            s_a[i] = i;
        }
        fclose(fp);
        return;
    }

    int opt;
    const char* optstring = "m:s:d:";

    while ((opt = getopt(argc, argv, optstring)) != -1) {
        FILE* fp = fopen(optarg, "r");
        if (opt == 'm') {
            for (int i = 0; i < n; i++) {
                fscanf(fp, "%lf%d", &a[i].v, &a[i].t);
                s_a[i] = i;
            }
        }
        fclose(fp);
    }
}

// 計算退火解的盈利
double calculate(SAnode s) {
    Koutput str;
    double profit = 0;
    ksj1.clear(), ksj2.clear();
    for (int i = 0; i < s.n; i++) {
        int u = s.B[i];
        if (a[u].t <= 10) str.d[str.n][0] = 0, str.d[str.n][1] = u, str.n++, profit -= a[u].v;
        else ksj2.pb(u);
    }
    // 初始化
    for (int i = 1; i <= 10; i++) {
        for (int j = 0; j < 3; j++) {
            if (a[s.A[i][j]].t <= 5) ksj1.pb(s.A[i][j]);  //=
            else ksj2.pb(s.A[i][j]);
        }
    }
    vector<int>::iterator it;
    // 開始遍歷每一天
    for (int i = 1; i <= 10; i++) {
        sort(s.A[i], s.A[i] + 3, cmp_greater);
        double maxv = a[s.A[i][0]].v;

        // 把臨近過期的藥轉移到 ksj1
        for (it = ksj2.begin(); it != ksj2.end(); ) {
            int u = *it;
            if (a[u].t - i + 1 <= 5) ksj1.pb(u), it = ksj2.erase(it);  //=
            else it++;
        }

        // 上架賣的三種藥
        for (int j = 0; j < 3; j++) {
            int u = s.A[i][j];
            str.s[i][j][0] = u;
            if (a[u].t - i + 1 <= 5) ksj1.erase(find(ksj1.begin(), ksj1.end(), u));  //=
            else ksj2.erase(find(ksj2.begin(), ksj2.end(), u));
        }
        sort(ksj1.begin(), ksj1.end(), cmp_greater);
        sort(ksj2.begin(), ksj2.end(), cmp_greater);

        profit -= ksj1.size() + ksj2.size() * 0.5;  // 先減去所有藥的管理費，這樣可以把上架藥品的管理費看做收益

        int p = 0, q = 0;
        int k0 = 6, k1 = 6, k2 = 6, u = s.A[i][0], v = s.A[i][1], w = s.A[i][2];
        int maxk0, maxk1, maxk2, maxp, maxq;
        double Max = -10000;

        // 列舉賣出藥品的定價不大於 maxv+k，計算此時當天收益（賣藥利潤+上架不賣的7件藥品的管理費）
        for (int k = 6; k >= -1; k--) {
            while (p < ksj1.size() && a[ksj1[p]].v + 6 > maxv + k) p++;
            while (q < ksj2.size() && a[ksj2[q]].v + 6 > maxv + k) q++;
            while (a[u].v + deltav[k0] > maxv + k) k0--;
            while (a[v].v + deltav[k1] > maxv + k) k1--;
            while (a[w].v + deltav[k2] > maxv + k) k2--;

            double tmp = deltav[k0] + deltav[k1] + deltav[k2] + min(p, 7) + 0.5 * min(q, 7 - min(p, 7));
            if (tmp > Max) {
                Max = tmp;
                maxp = p, maxq = q;
                maxk0 = k0, maxk1 = k1, maxk2 = k2;
            }
        }
        int t = 0;
        str.s[i][0][1] = maxk0;
        str.s[i][1][1] = maxk1;
        str.s[i][2][1] = maxk2;

        while (t < 7 && maxp) str.s[i][3 + t][0] = ksj1[--maxp], str.s[i][3 + t][1] = 6, t++;
        while (t < 7 && maxq) str.s[i][3 + t][0] = ksj2[--maxq], str.s[i][3 + t][1] = 6, t++;
        while (t < 7) str.s[i][3 + t][0] = -1, str.s[i][3 + t][1] = 6, t++;
        profit += Max;
    }
    
    if (profit > bestans) {
        bestans = profit;
        beststr = str;
        bestSAstr = s;
    }

    return profit;
}

// 調整退火解
SAnode adjust(SAnode s) {
    for (int i = 0; i <= rand() % 5; i++) {
        while (1) {
            int opt = rand() % 4;
            if (opt == 0) {
                int x = rand() % 10 + 1, y = rand() % 3, z = rand() % s.n;
                if (a[s.B[z]].t >= x && a[s.A[x][y]].t > 10) {
                    swap(s.B[z], s.A[x][y]);
                    break;
                }
            } else {
                int x = rand() % 10 + 1, y = rand() % 3;
                int u = rand() % 10 + 1, v = rand() % 3;
                if (a[s.A[x][y]].t >= u && a[s.A[u][v]].t >= x) {
                    swap(s.A[x][y], s.A[u][v]);
                    break;
                }
            }
        }
    }
    return s;
}

// 貪心生成初始解
SAnode greedy() {
    int p = 0;
    SAnode ret;
    for (int i = 1; i <= 10; i++) {
        for (int j = 0; j < 3; j++) {
            while (a[s_a[p]].t < i) {
                int u = s_a[p];
                for (int ii = 1; ii <= a[u].t; ii++) {
                    for (int jj = 0; jj < 3; jj++) {
                        if (a[ret.A[ii][jj]].v < a[u].v) swap(u, ret.A[ii][jj]);
                    }
                }
                initDel[initDelcnt++] = u;
                initCost += a[u].v;
                p++;
            }
            ret.A[i][j] = s_a[p++];
        }
    }
    while (p < 50) {
        ret.B[ret.n++] = s_a[p++];
    }
    return ret;
}

// 進行模擬退火
double SA() {
    double T = 10000, alpha = 0.99, lst = -10000;
    SAnode state = initstate;
    lst = calculate(state);
    while (T > 0.000001) {
        SAnode newState = adjust(state);
        double tmp = calculate(newState), delta = tmp - lst;
        if (delta > 0 || exp(delta/T) > frand()) lst = tmp, state = newState;
        T *= alpha;
    }
    return lst;
}

int main(int argc, char* argv[]) {
    srand(time(0) + 9927);
    //cout<<"#";
    read(argc, argv);

    sort(s_a, s_a + n, cmp_id);
    initstate = greedy();
    for (int i = 1; i <= 300; i++) SA();

    // 在所得解附近爬山
    SAnode state = bestSAstr;
    double lst = bestans;
    for (int i = 1; i <= 5000; i++) {
        SAnode newState = adjust(state);
        double tmp = calculate(newState);
        if (tmp > lst) {
            lst = tmp;
            state = newState;
        }
    }
    
    for (int i = 0; i < initDelcnt; i++) beststr.d[beststr.n][0] = 0, beststr.d[beststr.n][1] = initDel[i], beststr.n++; 
    beststr.print();
    cout << bestans - initCost;
    return 0;
}
/* by Hellsegamosken */