單純形法


看了集訓隊答辯,感覺要學習的有杜教篩高級版、線性規划、FFT、仙人掌、高級版線段樹

不出意外的話一個月內博客內都不會有別的東西了QAQ

首先是喜聞樂見的單純形法解線性規划。

今年(2016年)和線性規划有關的集訓隊論文有兩篇,大家可以自行翻一下集訓隊論文(當然如果你沒有拿到你可以去UOJ群下載啊),下面的大部分內容都是參閱akf那篇

線性規划的標准型一般長得像這樣:

image

一般我們拿到的都是像標准型這樣的問題,例如網絡流問題,我們就是要最大化流量並且讓每個點流量守恆。

假設我們把匯點到源點連一條容量為+∞的邊那么所有點就都流量守恆了。我們用c(u,v)表示(u,v)這條邊的容量,f(u,v)表示這條邊的流量。

那么我們不難得到這樣的一個標准型線性規划

image

不難把它轉化成標准型這樣的模型。

而最小費用流則也是一個標准型,我們還是把匯點到源點加一條容量為+∞,費用為0的邊。

那么也可以得到一個比較標准的線性規划:

image

(較論文中的配文有修改)

如果你要最小費用最大流的話只要加一個約束條件$\sum_{(u,v)}{f(u,v)}=maxflow$image

額那不等號怎么辦呢?

對於不等約束條件${\sum_{j=1}^n{a_{ij}x_j}}\leq{b_i}$image

在松弛型等式左邊的那些變量我們把它們叫做基變量,在上面的松弛型表示中就是x[n+1]…x[n+m],非基變量就是在等式右邊的那些,在上面的表示中是x[1]…x[n]。

顯然當bi非負時令所有基變量為bi,非基變量為0,即可得到一個滿足條件的初始解。(至於bi可能為負的情況下面在單獨說)

單純形法有一個重要的操作,叫轉軸(pivot)操作。就是說我們可以把一個基變量x[b]和非基變量x[n]互換,用x[b]和其他非基變量代換這個x[n],這樣x[b]就成了非基變量,而x[n]成了基變量。

一開始我們知道image

那么簡單代換一下會發現image,然后我們也可以把其他約束條件中的x[n]這樣代換掉,於是就完成了這個轉軸操作。

然后有了這個轉軸過程的幫助,我們就可以實現線性規划算法了。

首先為了最大化目標函數,我們考慮不停地找一個系數為正的非基變量,然后增大這個x。

具體的這個最優化過程如下:

image

(注意鄒逍遙原論文在這里的做法是選擇滿足image的,ysy告訴我這樣是不科學的,akf的論文里寫的也是ce>0,於是做了一些修改)

這是為什么呢?akf在論文中給出了一個例子,從中我們可以大概理解一下。

這個線性規划是這樣的:

image

我們考慮第一步我們選擇換正系數的x1,因為增大x1必然會增大目標函數。

我們分別考慮三個限制,顯然第一個限制下x1<=30,第二個限制下x1<=12,第三個限制x1<=9,然后x1<=9的下界最緊,所以我們用x1代換x6,得到下面的新線性規划:

image

哇目標函數增大到了27,可喜可賀。

接下來假設我們來增大x3,類似地可以得到:

image

然后可以增大x2,可以得到:

image

由於所有非基變量系數均為負的,所以我們不可能再得到更小的解了。

image

咦,有沒有注意到上面漏了什么東西說要下面講啊…

對了,在bi為負的時候,如果你把所有基變量帶0,而非基變量就不一定是一組合法的初始解。

這時候我們就需要一個初始化操作,初始化操作的基本思想是引入一個輔助線性規划:

image

如果我們求得了這個輔助線性規划的最優解,那么如果最優解中x0>0顯然原線性規划無解。

如果最優解中x0=0,那么x[1]…x[n+m]就是原線性規划的一個可行解。

而我們容易構造出輔助線性規划的一個可行解。

我們只要把x0作為換入變量,找到bi的最小值bl,把x[i+l]用x0來替代。

例如:

image

我們引入輔助線性規划:

image

bi最小值是-4,那么我們把x4用x0來替代:

image

然后我們就可以得到一組x的初始解用來進行最優化操作。

講道理:為什么bi這樣之后一定為正?

首先bi的最小值bl一定為負(廢話,否則就不用進行最優化操作了)

然后我們考慮第l個約束會變為:

image

而-bi顯然為正的。

而其他約束會變為:

image

由於bl是最小值,所以bi>=bl,所以-bl+bi>=0。

而我們發現UOJ上似乎有一種更加神奇的初始化方法?

我們的目標顯然是讓所有系數bi都為非負的。

我們選擇一個bi為負的基變量x[i+n],然后我們在該約束右邊找一個系數為正的非基變量,然后進行轉軸操作,如果沒有系數為正顯然就無解了。

額其實這和這個初始化操作本質上是一樣的。

所以這樣就可以完成整個線性規划的過程。

我們來分析一下時間復雜度?

pivot操作的復雜度顯然是O(NM)的,但是最優化操作中pivot操作的調用次數可能會成為指數級。

但是我們可以發現要達到這個指數級的調用次數,邊權也必須是指數級的。所以在OI中往往跑得比誰都快。所以“能在1s內跑出范圍為幾百的數據”。

然而一般線性規划是可以在多項式范圍內求解的,只不過…

橢球算法 O(n^6*m^2)

內點算法 O(n^3.5*m^2)

改進的內點算法 O(n^3.5*m)

在oi中一般並不能有什么卵用

所以單純形法就這么講完啦。

實現的時候可以發現,寫單純形並不要真的把松弛型建出來,只要假裝第i個約束條件就對應第i個基變量,把系數取負就可以了。

科學的代碼 http://uoj.ac/submission/61872 (worse)

我的代碼(有一些細節寫的比較丑

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <algorithm>
#include <string.h>
#include <vector>
#include <math.h>
#include <limits>
#include <set>
#include <map>
using namespace std;
typedef double ld;
#define SZ 233
int n,m,t,id[SZ];
ld a[SZ][SZ],vv[SZ];
const ld eps=1e-7;
int dcmp(ld x)
{
    if(x<-eps) return -1;
    if(x>eps) return 1;
    return 0;
}
void pivot(int r,int c) //基變量r和非基變量c 
{
    swap(id[r+n],id[c]);
    ld x=-a[r][c]; a[r][c]=-1;
    for(int i=0;i<=n;i++) a[r][i]/=x;
    for(int i=0;i<=m;i++)
    {
        if(dcmp(a[i][c])&&i!=r);else continue;
        x=a[i][c]; a[i][c]=0;
        for(int j=0;j<=n;j++) a[i][j]+=x*a[r][j];
    }
}
void solve()
{
    for(int i=1;i<=n;i++) id[i]=i;
    while(1) //init
    {
        int x=0,y=0;
        for(int i=1;i<=m;i++)
        {
            if(dcmp(a[i][0])<0&&(!x||(rand()&1))) x=i; 
        }
        if(!x) break;
        for(int i=1;i<=n;i++)
        {
            if(dcmp(a[x][i])>0&&(!y||(rand()&1))) {y=i; break;}
        }
        if(!y) {puts("Infeasible"); return;}
        pivot(x,y);
    }
    while(1) //simplex
    {
        int x=0,y=0;
        for(int i=1;i<=n;i++)
        {
            if(dcmp(a[0][i])>0&&(!x||(rand()&1))) x=i; 
        }
        if(!x) break;
        double w,t; bool f=1;
        for(int i=1;i<=m;i++)
        {
            if(dcmp(a[i][x])<0&&((t=-a[i][0]/a[i][x]),f||t<w)) w=t, y=i, f=0;
        }
        if(!y) {puts("Unbounded"); return;}
        pivot(y,x);
    }
    printf("%.9lf\n",a[0][0]);
    if(!t) return;
    for(int i=1;i<=n;i++) vv[i]=0;
    for(int i=n+1;i<=n+m;i++) vv[id[i]]=a[i-n][0];
    for(int i=1;i<=n;i++) printf("%.9lf ",vv[i]);
}
int main()
{
    scanf("%d%d%d",&n,&m,&t);
    for(int i=1;i<=n;i++) scanf("%lf",&a[0][i]);
    for(int i=1;i<=m;i++)
    {
        for(int j=1;j<=n;j++) scanf("%lf",&a[i][j]), a[i][j]*=-1;
        scanf("%lf",&a[i][0]);
    }
    solve();
}


免責聲明!

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



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM