uoj#179 線性規劃 單純形法の模板
這是一道模板題。
本題中你需要求解一個標準型線性規劃:
有n個實數變數x1,x2,⋯,xn和m條約束,其中第i條約束形如aij*xj≤bi ,j∈(1,n),i∈(1,m)
此外這n個變數需要滿足非負性限制,即xj≥0。
在滿足上述所有條件的情況下,你需要指定每個變數xj的取值,使得目標函式F=cj*xj ,j∈(1,n)的值最大。輸入格式
第一行三個正整數 n,m,t。其中t∈{0,1}。
第二行有n個整數c1,c2,⋯,cn,整數間均用一個空格分隔。
接下來m行,每行代表一條約束,其中第i行有n+1個整數ai1,ai2,⋯,ain,bi,整數間均用一個空格分隔。
輸出格式
如果不存在滿足所有約束的解,僅輸出一行”Infeasible”。
如果對於任意的M,都存在一組解使得目標函式的值大於M,僅輸出一行”Unbounded”。
否則,第一行輸出一個實數,表示目標函式的最大值F。當第一行與標準答案的相對誤差或絕對誤差不超過10−6,你的答案被判為正確。
如果t=1,那麼你還需要輸出第二行,用空格隔開的n個非負實數,表示此時x1,x2,⋯,xn的取值,如有多組方案請任意輸出其中一個。
判斷第二行是否合法時,我們首先檢驗F−cjxj,j∈(1,n)是否為0,再對於所有ii,檢驗min{0,bi−aijxj,j∈(1,n)}是否為0。檢驗時我們會將其中大於0的項和不大於0的項的絕對值分別相加得到S+和S−,如果S+和S−的相對誤差或絕對誤差不超過10−6,則判為正確。
如果t=0,或者出現Infeasible或Unbounded時,不需要輸出第二行。
-update 2017/3/20
對了,對於是否有初始可行解的話,initialize這個函式是一種方法,但是可能會有點慢。大多時候都利用對偶原理來轉化模型。額奧爺爺說當c[]都是負的時候才用對偶(?)
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<ctime>
#include<iostream>
#include<algorithm>
using namespace std;
#define maxn 25
const double eps=0.00000001 ,inf=1e15;
int n,m,id[maxn*2];
double a[maxn][maxn];//a[i][j]:i表第幾條約束 j表第幾個元素
double myabs(double x) {return x>0?x:-x;}
//a[0][i] -> ci 目標函式中第i個元素係數
//a[i][0] -> bi 第i條約束中的常數
//a[i][j] -> Aij 第i條約束中第j個元素的係數
//最大化 sigma(ci*xi),i∈N
//約束 xj=bj-sigma(aji*xi) ,j∈B
//轉軸
void pivot(int l,int e)
//替入變數xe∈非基本變數(1~n) 替出變數xl∈基本變數(n+1~n+m)
{
int tt=id[n+l];id[n+l]=id[e];id[e]=tt;
int i,j;double t=a[l][e];a[l][e]=1;
for (j=0;j<=n;j++) a[l][j]/=t;
//重寫其它等式:
for (i=0;i<=m;i++)
if (i!=l && myabs(a[i][e])>eps)
{
t=a[i][e];a[i][e]=0;
for (j=0;j<=n;j++)
a[i][j]-=a[l][j]*t;
}
}
//初始化
//方法一:引入一個輔助線性規劃 要求最大化-x0
//約束為 xj=bj-sigma(aji*xi)+x0 ,j∈B然後用x0替換bj為負的約束
//下面的是方法二:
bool initialize()
{
while (1)
{
int i,j,e=0,l=0;
for (i=1;i<=m;i++)
if (a[i][0]<-eps && (!l || (rand()&1))) l=i;
if (!l) break;
for (j=1;j<=n;j++)
if (a[l][j]<-eps && (!e || (rand()&1))) e=j;
if (!e) {printf("Infeasible\n");return 0;}
pivot(l,e);
//在bi為負的時候,把所有基變數設為0不是一組合法的初始解
//所以選擇一個bi為負的基變數x[i+n]
//然後在該約束右邊找一個係數為正(即原係數為負)的非基變數進行轉軸操作
//如果沒有係數為正顯然就無解了
}return 1;
}
//最優化
bool simplex()
{
int i,j;
while (1)
{
int l=0,e=0;double minn=inf;
for (j=1;j<=n;j++)
if (a[0][j]>eps) {e=j;break;}
if (!e) break;
//如果目標變數ci都小於0 那麼最優解就是非基變數都為0
for (i=1;i<=m;i++)
if (a[i][e]>eps && a[i][0]/a[i][e]<minn)
minn=a[i][0]/a[i][e],l=i;
//在所有的式子中找出包含當前選中項(係數不為0)且最緊的一項
if (!l) {printf("Unbounded\n");return 0;}
//如果所有的a[i][e]都小於0,說明最優值正無窮
pivot(l,e);
}return 1;
}
double ans[maxn];
int main()
{
//freopen("a.in","r",stdin);
//freopen("a.out","w",stdout);
srand(time(0));int t,i,j;
scanf("%d%d%d",&n,&m,&t);
//n個變數 m條約束
for (i=1;i<=n;i++) scanf("%lf",&a[0][i]);
for (i=1;i<=m;i++)
{
for (j=1;j<=n;j++)
scanf("%lf",&a[i][j]);
scanf("%lf",&a[i][0]);
}
for (i=1;i<=n;i++) id[i]=i;
if (initialize() && simplex())
{
printf("%.8lf\n",-a[0][0]);
if (t)
{
for (i=1;i<=m;i++) ans[id[n+i]]=a[i][0];
for (i=1;i<=n;i++) printf("%.8lf ",ans[i]);
}
}
return 0;
}