【BZOJ3265】志願者招募加強版 線性規劃 單純形法 對偶原理

空灰冰魂發表於2015-05-12

連結:

#include <stdio.h>
int main()
{
    puts("轉載請註明出處[vmurder]謝謝");
    puts("網址:blog.csdn.net/vmurder/article/details/45673079");
}

題解:

這道題是線性規劃求目標函式最小值 ,對偶原理轉一下就成了單純形演算法求線性規劃最大值。

單純形法:

首先這篇部落格缺失了很多證明,只能講述單純形法的實現。

 // N個變數 M個限制
double a[M][N]; // 這m個限制裡變數的係數
double b[M]; // 第i個限制加和 <=bi
double c[N]; // 目標函式裡變數係數
double ans; // 目標函式自帶的係數
void pivot(int x,int y); // 轉軸操作
void solve();

首先對於如下資料

3 3 

2 3 4

1 1 2 2

1 2 3 5

1 3 3 2

我們可以得到以下的限制:
目標函式: 2x1+5x2+2x3+0

2x_1+5x_2+2x_3+0

限制函式:
x12
x_1\leq2

x1+x23
x_1+x_2\leq3

x2+x34
x_2+x_3\leq4

然後若 x4x5x6

x_4 x_5 x_6
均為任意值。
那麼我們可以將式子表示為
狀態 1: **************
目標函式: 2x1+5x2+2x3+0
2x_1+5x_2+2x_3+0

限制函式:
x4+1x1=2
x_4+1x_1=2

x5+1x1+1x2=3
x_5+1x_1+1x_2=3

x6+1x2+1x3=4
x_6+1x_2+1x_3=4

然後在單純形的過程中函式們的各系數得到改變的過程如下:

狀態 1: **************
目標函式: 2x1+5x2+2x3+0

2x_1+5x_2+2x_3+0

限制函式:
x4+1x1=2
x_4+1x_1=2

x5+1x1+1x2=3
x_5+1x_1+1x_2=3

x6+1x2+1x3=4
x_6+1x_2+1x_3=4

狀態 2: **************
目標函式: 2x4+5x2+2x3+4
-2x_4+5x_2+2x_3+4

限制函式:
x1+1x4=2
x_1+1x_4=2

x5+1x4+1x2=1
x_5+-1x_4+1x_2=1

x6+1x2+1x3=4
x_6+1x_2+1x_3=4

狀態 3: **************
目標函式: 3x4+5x5+2x3+9
3x_4+-5x_5+2x_3+9

限制函式:
x1+1x4=2
x_1+1x_4=2

x2+1x4+1x5=1
x_2+-1x_4+1x_5=1

x6+1x4+1x5+1x3=3
x_6+1x_4+-1x_5+1x_3=3

狀態 4: **************
目標函式: 3x1+5x5+2x3+15
-3x_1+-5x_5+2x_3+15

限制函式:
x4+1x1=2
x_4+1x_1=2

x2+1x1+1x5=3
x_2+1x_1+1x_5=3

x6+1x1+1x5+1x3=1
x_6+-1x_1+-1x_5+1x_3=1

狀態 5: **************
目標函式: 1x1+3x5+2x6+17
-1x_1+-3x_5+-2x_6+17

限制函式:
x4+1x1=2
x_4+1x_1=2

x2+1x1+1x5=3
x_2+1x_1+1x_5=3

x3+1x1+1x5+1x6=1
x_3+-1x_1+-1x_5+1x_6=1

線性規劃最終答案:17

單純形思想:

其實就是每次找一個在目標函式中存在的變數,跟不在其中的某個變數交換一下,然後換來換去最終目標函式中所有變數係數都為負了,那麼目標函式最後加的常數就是答案辣。

單純形實現:

具體的實現過程呢?

struct Simplex
{
    int n,m; // n個變數、m個限制
    double a[M][N],b[M],c[M],ans;
    void pivot(int x,int y)
    {
        int i,j;
        double t;
        b[x]/=a[x][y];
        for(i=1;i<=n;i++)if(i!=y)a[x][i]/=a[x][y];
        a[x][y]=1.0/a[x][y];

        for(i=1;i<=m;i++)if(i!=x&&fabs(a[i][y])>eps)
        {
            b[i]-=a[i][y]*b[x],t=a[i][y];
            for(j=1;j<=n;j++)a[i][j]-=t*a[x][j];
            a[i][y]=-t*a[x][y];
        }

        ans+=c[y]*b[x],t=c[y];
        for(i=1;i<=n;i++)c[i]-=t*a[x][i];
        c[y]=-t*a[x][y];
    }
    double solve()
    {
        read();
        double t;
        for(int i,x,y;;)
        {
            for(i=1;i<=m;i++)if(c[i]>eps){y=i;break;}
            if(i>m)return ans;

            for(t=inf,i=1;i<=m;i++)
                if(a[i][y]>eps&&t>b[i]/a[i][y])
                    x=i,t=b[i]/a[i][y];
            if(t==inf)return t;
            else pivot(x,y);
        }
    }
}simplex;

下面的實現裡存在兩個坑,請不要深究它們怎麼證。

double solve()

double~solve()
:

首先我們每次找一個標號最靠前的【坑1】目標函式中係數為正的變數 y

y
(如第30行),然後如果找不到,就表示已經找到了最優解【坑2】。

然後我們找這個變數在哪個限制裡是”最緊”的【坑3】,即33、34、35這三行。
如果找不到任何一個”緊”的行,即 a(i,y)<=eps

\forall a(i,y)<=eps
那麼這個變數想多大就能多大,反正有輔助變數頂著,可以直接返回無界【坑4】,即答案是無窮大。

然後我們對此限制中的輔助變數和找到的 y

y
進行轉軸操作。

void pivot(int x,int y)

void~pivot(int~x,int~y)
:

先把第 x

x
個限制的係數改一改,然後用它去消其它的限制,再消一下目標函式。結束。類似高斯消元?

對偶原理:

如果要求最小值,那麼我們把模型對偶一下就可以辣。【坑5】

即限制矩陣 【轉置】 一下,然後目標函式的係數和每個限制的最終

\leq
的那個值也互換一下。

【轉置】 :矩陣的元素 a(i,j)

a(i,j)
變為 a(j,i)
a(j,i)

對於證明坑的總結:

【坑1】為什麼要找標號最小的?
根據Bland法則,這麼做可以避免被卡死迴圈。
【坑2】為什麼標號都是負的就表示找到了最優解?
表示無法理解為什麼不能給某個存在負係數的變數轉一轉。
【坑3】為什麼要找”最緊”的?
不知道不知道。
【坑4】為什麼直接返回無界了?
不可以別的某個變數再來轉一轉,然後就迴歸有界?
【坑5】對偶原理是毛線?
天然坑,深坑。

程式碼:

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define N 1010 // 變數個數
#define M 10100 // 限制個數
#define inf 1e9
#define eps 1e-5
using namespace std;
struct Simplex
{
    int n,m; // n個變數、m個限制
    double a[M][N],b[M],c[M],ans;
    void read()
    {
        int i,j,k,l,r;
        scanf("%d%d",&m,&n);
        for(i=1;i<=m;i++)scanf("%lf",&c[i]);
        for(i=1;i<=n;i++)
        {
            scanf("%d",&k);
            while(k--)
            {
                scanf("%d%d",&l,&r);
                for(j=l;j<=r;j++)a[i][j]=1.0;
            }
            scanf("%lf",&b[i]);
        }
        swap(n,m);
    }
    void pivot(int x,int y)
    {
        int i,j;
        double t;
        b[x]/=a[x][y];
        for(i=1;i<=n;i++)if(i!=y)a[x][i]/=a[x][y];
        a[x][y]=1.0/a[x][y];

        for(i=1;i<=m;i++)if(i!=x&&fabs(a[i][y])>eps)
        {
            b[i]-=a[i][y]*b[x],t=a[i][y];
            for(j=1;j<=n;j++)a[i][j]-=t*a[x][j];
            a[i][y]=-t*a[x][y];
        }

        ans+=c[y]*b[x],t=c[y];
        for(i=1;i<=n;i++)c[i]-=t*a[x][i];
        c[y]=-t*a[x][y];
    }
    double solve()
    {
        read();
        double t;
        for(int i,x,y;;)
        {
            for(i=1;i<=m;i++)if(c[i]>eps){y=i;break;}
            if(i>m)return ans;

            for(t=inf,i=1;i<=m;i++)
                if(a[i][y]>eps&&t>b[i]/a[i][y])
                    x=i,t=b[i]/a[i][y];
            if(t==inf)return t;
            else pivot(x,y);
        }
    }
}simplex;
int main()
{
//  freopen("test.in","r",stdin);
    printf("%d\n",(int)(simplex.solve()+eps));
    return 0;
}

相關文章