最小费用最大流问题(看了一周《算法导论》的单纯性算法已经长跪不起了!):在最大流的基础上,使各个路径的花费和最小。
第一种方法:参考最大流。
简要分析一下两种问题的异同,最大流问题建图时只有有向边的容量capacity,而最小费用最大流还有有向边的单位容量的花费cost,所以找增广路径时略有不同。
- 简单地说,最大流的增广路径的可行边是有流量的边就可以走,能到终点就是增广路径。
- 而最小费用最大流的增广路径必须是:有容量且单位容量花费和最小的路径。
为什么是单位容量最小花费和,而不是增广路径最小总花费和(路径最小容量x路径单位容量花费和),其实易知第二个值是不容易求的,因为它不是最短路而是最小乘法加权值之和。具体原理待证明。
算法具体实现:对每件商品建图求最小费用最大流
- 建图:新建超级源节点到所有供应商的容量为供应商的存储量,同时供应商到商店的容量初值也设置为供应商存储量,新建超级汇点,商店到汇点容量为商店需求量;供应商到商店花费为单位商品运费,同时取其相反数赋给商店到供应商花费(相当于退货)。
- 找图中增广路径,若失败则说明已最优;反之根据路径更改容量(同最大流方法)。
求最小单位容量花费路径可以使用SPFA,Bellman_Ford。
第二种方法:
KM算法
KM算法是一种求最大加权二分匹配的方法,运用于此题,首先注意到二分匹配图中权值可以等效于单位容量花费,但是其中没有边的容量概念。但是我们可以把商店和供应商分解成若干(需求量或者供应量)相当于容量为1的小商店以及小供应商,就可以等价于求二分最大加权匹配了。
算法具体实现:
- 根据需求量把商店分解到二分图中X集合中的点,同理根据供应量把供应商分解为二分图中Y集合中的点,各边权值仍然为原商店与原供应商的单位花费的相反数。可知X中的点一定是有匹配的,即需求是必须满足的,但Y中的点不一定有匹配(供大于求)。
- 根据X中点序,直到能找到一条增广路径再继续X中下一个点,若找不到则修改左右标签值。
SPFA:
#include<iostream>
#include<cstring>
#include<cstdio>
#include<queue>
using namespace std;
#define INF 0X7F7F
#define maxN 55
#define maxK 55
int N,M,K;
short need[maxK][maxN];
short storage[maxK][maxN];
short capacity[(maxN)<<1][(maxN)<<1];
short cost[(maxN)<<1][(maxN)<<1];
short minCost;
short dest;
short spfa(short* fa)
{
short i,now;
short minDis[(maxN)<<1];
char vis[(maxN)<<1];
memset(vis,0,sizeof(vis));
memset(minDis,127,sizeof(minDis));
minDis[0] = 0;
queue<short> q;
q.push(0);
vis[0] = true;
while(!q.empty())
{
now = q.front();
q.pop();
vis[now] = false;
for(i = 1;i <= dest;i++)
{
if(capacity[now][i] > 0&&minDis[now]+cost[now][i] < minDis[i]) //找出有容量且单位商品的最小花费路径
{
fa[i] = now;
minDis[i] = minDis[now]+cost[now][i];
if(!vis[i])
{
q.push(i);
vis[i] = true;
}
}
}
}
return minDis[dest];
}
int minCost_maxFlow()
{
short x;
short minCapacity = INF;
short fa[(maxN)<<1];
memset(fa,0,sizeof(fa));
while(spfa(fa) != INF)
{
x = dest;
while(x) //找出增广路径的最小容量
{
minCapacity = min(minCapacity,capacity[fa[x]][x]);
x = fa[x];
}
x = dest;
while(x) //修改增广路径中容量
{
minCost += minCapacity*cost[fa[x]][x];
capacity[fa[x]][x] -= minCapacity;
capacity[x][fa[x]] += minCapacity;
x = fa[x];
}
memset(fa,0,sizeof(fa));
}
return 0;
}
int main()
{
int i,j,k;
short sumNeed[maxK];
short sumStorage[maxK];
char needEnough;
while(~scanf("%d%d%d",&N,&M,&K)&&N)
{
minCost = 0;
dest = N+M+1;
needEnough = true;
memset(cost,0,sizeof(cost));
memset(sumNeed,0,sizeof(sumNeed));
memset(sumStorage,0,sizeof(sumStorage));
for(i = 1;i <= N;i++)
{
for(k = 1;k <= K;k++)
{
scanf("%d",&need[k][i]); //第k个商家需要多少第j种商品
sumNeed[k] += need[k][i];
}
}
for(i = 1;i <= M;i++)
{
for(k = 1;k <= K;k++)
{
scanf("%d",&storage[k][i]); //第k个供货商有多少j的库存
sumStorage[k] += storage[k][i];
}
}
for(k = 1;k <= K;k++)
{
if(sumNeed[k] > sumStorage[k]) //第k中商品的总供应量少于需求
{
needEnough = false;
break;
}
}
for(k = 1;k <= K;k++)
{
memset(capacity,0,sizeof(capacity));
for(i = 1;i <= N;i++)
{
capacity[M+i]dest] = need[k][i]; //商店到汇点的容量为商店的需求
for(j = 1;j <= M;j++)
{
capacity[j][M+i] = capacity[0][j] = storage[k][j]; //源点到供应商以及供应商到商店的容量为供应商存储量
scanf("%d",&cost[j][M+i]);
cost[M+i][j] = -1*cost[j][M+i]; //商店到供应商的花费为供应商到商店花费的相反数(相当于退回商品)
}
}
if(needEnough)
minCost_maxFlow();
}
if(!needEnough)
{
printf("-1\n");
}
else
printf("%d\n",minCost);
}
return 0;
}
#include<iostream>
#include<cstring>
#include<cstdio>
#include<queue>
using namespace std;
#define INF 0X7F7F
#define maxN 55
#define maxK 55
#define maxE 5100
int N,M,K;
short need[maxK][maxN];
short storage[maxK][maxN];
short capacity[(maxN)<<1][(maxN)<<1];
short cost[(maxN)<<1][(maxN)<<1];
short edges[maxE][2];
short edgeNum;
short minCost;
short dest;
int generateEdge() //每次容量改变后都重新生成边集合
{
edgeNum = 0;
short i,j;
for(i = 0;i < dest;i++)
{
for(j = 1;j <= dest;j++)
{
if(capacity[i][j] > 0)
{
edges[edgeNum][0] = i;
edges[edgeNum++][1] = j;
}
}
}
return 1;
}
short bellman_ford(short* fa)
{
short i,j;
short x,y;
short minDis[(maxN)<<1];
char IsRelax;
memset(minDis,0x7F,sizeof(minDis));
minDis[0] = 0;
for(i = 0;i < dest;i++)
{
IsRelax = false;
for(j = 0;j < edgeNum;j++)
{
x = edges[j][0];
y = edges[j][1];
if(minDis[y] > minDis[x]+cost[x][y])
{
fa[y] = x;
minDis[y] = minDis[x]+cost[x][y];
IsRelax = true;
}
}
if(!IsRelax)
break;
}
return minDis[dest];
}
int minCost_maxFlow()
{
short x;
short minCapacity;
short fa[(maxN)<<1];
while(generateEdge()&&bellman_ford(fa) != INF)
{
minCapacity = INF;
x = dest;
while(x)
{
minCapacity = min(minCapacity,capacity[fa[x]][x]);
x = fa[x];
}
x = dest;
while(x)
{
minCost += minCapacity*cost[fa[x]][x];
capacity[fa[x]][x] -= minCapacity;
capacity[x][fa[x]] += minCapacity;
x = fa[x];
}
memset(fa,0,sizeof(fa));
}
return 0;
}
int main()
{
int i,j,k;
short sumNeed[maxK];
short sumStorage[maxK];
char needEnough;
while(~scanf("%d%d%d",&N,&M,&K)&&N)
{
minCost = 0;
dest = N+M+1;
needEnough = true;
memset(cost,0,sizeof(cost));
memset(sumNeed,0,sizeof(sumNeed));
memset(sumStorage,0,sizeof(sumStorage));
for(i = 1;i <= N;i++)
{
for(k = 1;k <= K;k++)
{
scanf("%d",&need[k][i]); //第k个商家需要多少第j种商品
sumNeed[k] += need[k][i];
}
}
for(i = 1;i <= M;i++)
{
for(k = 1;k <= K;k++)
{
scanf("%d",&storage[k][i]); //第k个供货商有多少j的库存
sumStorage[k] += storage[k][i];
}
}
for(k = 1;k <= K;k++)
{
if(sumNeed[k] > sumStorage[k])
{
needEnough = false;
break;
}
}
for(k = 1;k <= K;k++)
{
memset(capacity,0,sizeof(capacity));
for(i = 1;i <= N;i++)
{
capacity[M+i][M+N+1] = need[k][i];
for(j = 1;j <= M;j++)
{
capacity[j][M+i] = capacity[0][j] = storage[k][j];
scanf("%d",&cost[j][M+i]);
cost[M+i][j] = -1*cost[j][M+i];
}
}
if(needEnough)
minCost_maxFlow();
}
if(!needEnough)
{
printf("-1\n");
}
else
printf("%d\n",minCost);
}
return 0;
}
双向循环链表&Bellman_Ford(不用每次都生成边的集合):
#include<iostream>
#include<cstring>
#include<cstdio>
#include<queue>
using namespace std;
#define INF 0X7F7F
#define maxN 55
#define maxK 55
struct node
{
short nodeNo;
node* prev;
node* next;
};
node* adj[(maxN)<<1];
int N,M,K;
short need[maxK][maxN];
short storage[maxK][maxN];
short capacity[(maxN)<<1][(maxN)<<1];
short cost[(maxN)<<1][(maxN)<<1];
short minCost;
short dest;
int adjInsert(short a,short b)
{
node* newNode = new node;
newNode->nodeNo = b;
adj[a]->next->prev = newNode;
newNode->next = adj[a]->next;
newNode->prev = adj[a];
adj[a]->next = newNode;
return 0;
}
int adjInit()
{
short i,j;
for(i = 0;i < dest;i++) //初始化双向循环链表
{
adj[i] = new node;
adj[i]->prev = adj[i]->next = adj[i];
}
for(i = 0;i < dest;i++)
{
for(j = 1;j <= dest;j++)
{
if(capacity[i][j] > 0)
{
adjInsert(i,j);
}
}
}
return 0;
}
int adjDelete(short a,short b)
{
node* tmpNode = adj[a]->next;
while(tmpNode->nodeNo != b)
{
tmpNode = tmpNode->next;
}
tmpNode->prev->next = tmpNode->next;
tmpNode->next->prev = tmpNode->prev;
return 0;
}
short bellman_ford(short* fa)
{
short i,j;
short x;
node* tmpNode;
short minDis[(maxN)<<1];
char IsRelax;
memset(minDis,0x7F,sizeof(minDis));
minDis[0] = 0;
for(i = 0;i < dest;i++)
{
IsRelax = false;
for(j = 0;j < dest;j++)
{
for(tmpNode = adj[j]->next;tmpNode != adj[j];tmpNode = tmpNode->next)
{
x = tmpNode->nodeNo;
if(minDis[x] > minDis[j]+cost[j][x])
{
fa[x] = j;
minDis[x] = minDis[j]+cost[j][x];
IsRelax = true;
}
}
}
if(!IsRelax)
break;
}
return minDis[dest];
}
int minCost_maxFlow()
{
short x;
short minCapacity;
short fa[(maxN)<<1];
while(bellman_ford(fa) != INF)
{
minCapacity = INF;
x = dest;
while(x)
{
minCapacity = min(minCapacity,capacity[fa[x]][x]);
x = fa[x];
}
x = dest;
while(x)
{
minCost += minCapacity*cost[fa[x]][x];
if((capacity[fa[x]][x] -= minCapacity) == 0)
{
adjDelete(fa[x],x);
}
if(x != dest)
{
adjInsert(x,fa[x]);
}
capacity[x][fa[x]] += minCapacity;
x = fa[x];
}
}
return 0;
}
int main()
{
int i,j,k;
short sumNeed[maxK];
short sumStorage[maxK];
char needEnough;
while(~scanf("%d%d%d",&N,&M,&K)&&N)
{
minCost = 0;
dest = N+M+1;
needEnough = true;
memset(cost,0,sizeof(cost));
memset(sumNeed,0,sizeof(sumNeed));
memset(sumStorage,0,sizeof(sumStorage));
for(i = 1;i <= N;i++)
{
for(k = 1;k <= K;k++)
{
scanf("%d",&need[k][i]); //第k个商家需要多少第j种商品
sumNeed[k] += need[k][i];
}
}
for(i = 1;i <= M;i++)
{
for(k = 1;k <= K;k++)
{
scanf("%d",&storage[k][i]); //第k个供货商有多少j的库存
sumStorage[k] += storage[k][i];
}
}
for(k = 1;k <= K;k++)
{
if(sumNeed[k] > sumStorage[k])
{
needEnough = false;
break;
}
}
for(k = 1;k <= K;k++)
{
memset(capacity,0,sizeof(capacity));
for(i = 1;i <= N;i++)
{
capacity[M+i][M+N+1] = need[k][i];
for(j = 1;j <= M;j++)
{
capacity[j][M+i] = capacity[0][j] = storage[k][j];
scanf("%d",&cost[j][M+i]);
cost[M+i][j] = -1*cost[j][M+i];
}
}
if(needEnough)
{
adjInit();
minCost_maxFlow();
}
}
if(!needEnough)
{
printf("-1\n");
}
else
printf("%d\n",minCost);
}
return 0;
}
KM算法:
#include<iostream>
#include<cstring>
#include<cstdio>
using namespace std;
#define INF 0X7F7F
#define maxN 55
#define maxK 55
#define maxM 155
short cost[maxM][maxM]; //二分图左右权值
short lx[maxM],ly[maxM]; //二分图左右标签
char visx[maxM],visy[maxM]; //增广路中访问情况
short match[maxM]; //二分图匹配情况
short slack[maxM]; //松弛便签
short nx,ny;
short minCost;
bool DFS(short x)
{
short i,diff;
visx[x] = true;
for(i = 1;i <= ny;i++)
{
if(!visy[i])
{
diff = lx[x]+ly[i]-cost[x][i];
if(diff == 0) //增广路中可行边
{
visy[i] = true;
if(match[i] == 0||DFS(match[i])) //直到最后遇到Y中一点没有匹配才返回为1
{
match[i] = x;
return 1;
}
}
else if(slack[i] > diff) //i不属于增广路中Y的集合,更改其便签松弛量
{
slack[i] = diff;
}
}
}
return 0;
}
int KM()
{
short i,j;
short d;
memset(match,0,sizeof(match));
memset(lx,0X80,sizeof(lx));
memset(ly,0,sizeof(ly));
for(i = 1;i <= nx;i++) //初始化标签满足权值要求
{
for(j = 1;j <= ny;j++)
{
if(lx[i] < cost[i][j])
{
lx[i] = cost[i][j];
}
}
}
for(i = 1;i <= nx;i++)
{
memset(slack,0X7F,sizeof(slack)); //初始化松弛量为极大值
while(1)
{
memset(visx,0,sizeof(visx));
memset(visy,0,sizeof(visy));
if(DFS(i)) //找到增广路径就跳出循环,没找到改标号
{
break;
}
d = INF;
for(j = 1;j <= ny;j++) //找出使Y集合不在增广路中的点能满足要求的最大可松弛量d
{
if(!visy[j]&&slack[j] < d)
{
d = slack[j];
}
}
for(j = 1;j <= nx;j++) //X集合中在增广路中的点减去d
{
if(visx[j])
{
lx[j] -= d;
}
}
for(j = 1;j <= ny;j++)
{
if(visy[j]) //Y中集合在增广路中点加上d
{
ly[j] += d;
}
else //X中集合在增广路的点减去d,则Y中不在增广路中点松弛量还可减去d
{
slack[j] -= d;
}
}
}
}
for(i = 1;i <= ny;i++)
{
if(match[i]) //可能供大于求,Y中的点不一定有匹配
{
minCost += cost[match[i]][i];
}
}
return 0;
}
int main()
{
int N,M,K;
int i,j,k;
int tmpCost[maxN][maxN];
int need[maxK][maxN];
int storage[maxK][maxN];
short node2shop[maxM],node2supply[maxM];
short sumNeed[maxK];
short sumStorage[maxK];
char needEnough;
while(~scanf("%d%d%d",&N,&M,&K)&&N)
{
minCost = 0;
needEnough = true;
memset(sumNeed,0,sizeof(sumNeed));
memset(sumStorage,0,sizeof(sumStorage));
for(i = 1;i <= N;i++)
{
for(k = 1;k <= K;k++)
{
scanf("%d",&need[k][i]); //第k个商家需要多少第j种商品
sumNeed[k] += need[k][i];
}
}
for(i = 1;i <= M;i++)
{
for(k = 1;k <= K;k++)
{
scanf("%d",&storage[k][i]); //第k个供货商有多少j的库存
sumStorage[k] += storage[k][i];
}
}
for(k = 1;k <= K;k++)
{
if(sumNeed[k] > sumStorage[k])
{
needEnough = false;
break;
}
}
for(k = 1;k <= K;k++)
{
for(i = 1;i <= N;i++)
{
for(j = 1;j <= M;j++)
{
scanf("%d",&tmpCost[i][j]);
}
}
if(needEnough)
{
nx = ny = 1;
for(i = 1;i <= N;i++)
{
for(j = 1;j <= need[k][i];j++) //按照需求量将商店分解成若干需求为1的商店
{
node2shop[nx++] = i; //记录新生成的商店原来属于哪个商店
}
}
for(i = 1;i <= M;i++)
{
for(j = 1;j <= storage[k][i];j++) //按照供应量将供应商分解成若干供应量为1的供应商
{
node2supply[ny++] = i; //记录新生成的供应商原来属于哪个供应商
}
}
nx--;
ny--;
for(i = 1;i <= nx;i++)
{
for(j = 1;j <= ny;j++)
{
cost[i][j] = -1*tmpCost[node2shop[i]][node2supply[j]]; //转换成相反数求最大二分加权匹配
}
}
KM();
}
}
if(!needEnough)
{
printf("-1\n");
}
else
{
printf("%d\n",-1*minCost);
}
}
return 0;
}