对偶线性规划:n维参数向量满足两个等式,求另一个参数向量的极大值?
3.6与平面和空间打交道的计算几何
凸包
首先形式化描述该线性规划(线性规划与对偶问题的一般化详见附录):
原始问题:
max Z = C1X1 + C2X2 + C3X3 + … + CnXn
st A1X1 + A2X2 + A3X3 + … + AnXn = Si
B1X1 + B2X2 + B3X3 + … + BnXn = Ti
Xi ≥ 0
等效为:
max Z = C1X1 + C2X2 + C3X3 + … + CnXn
st A1X1 + A2X2 + A3X3 + … + AnXn ≤ Si
-A1X1 – A2X2 – A3X3 – … – AnXn ≤ –Si
B1X1 + B2X2 + B3X3 + … + BnXn ≤ Ti
-B1X1 – B2X2 – B3X3 – … – BnXn ≤ –Ti
Xi ≥ 0
但愿没人问出“将一个不等式左右两边乘以-1的话,不等符号应该反过来,第2个不等式为什么不是≥”这种问题,这么做是因为,a≥b && -a≥-b等效于a≥b && a≤b等效于a==b。
对偶形式:
min G = SiW1 – SiW2 + TiW3 – TiW4
st A1W1 – A1W2 + B1W3 – B1W4 ≥ C1
A2W1 – A2W2 + B2W3 – B2W4 ≥ C2
…
AnW1 – AnW2 + BnW3 – BnW4 ≥ Cn
Wi ≥ 0
令X=W1 – W2 ,Y=W3 – W4,则对偶形式简化为:
min G = SiX + TiY
st A1X + B1Y ≥ C1
A2X + B2Y ≥ C2
…
AnX + BnY ≥ Cn
于是我们得到n条直线,每一条都是一个约束。每一条直线将平面分割为两部分,只有其中一个半平面才是合乎约束的。这n个半平面的交集构成解空间,剩下的就是初中水平的一元二次线性规划了。
我是很想这么说,但是写起来真的好麻烦。
首先,一个严峻的问题是,半平面交集算法。
看到一种算法是,先按直线的方向向量角度和截距排序,然后维护一个栈,先压入前两条线,从第3条线开始,判断这条线与栈顶下一条线的交点是否满足栈顶的线,如果满足,弹出栈顶,继续判断,直到栈中元素不足两个,插入当前直线,继续加入下一条直线。画个图示意一下:
如图,直线标号表示加入顺序,红色表示弹出,绿色表示不弹出。
然后半平面交集有了后,对输入直线集合也排个序,在交集第一条直线左边的肯定无解,在最后一条直线右边的也无解。夹在中间的话,找角度最相似(最接近平行,也就是叉积最接近0)的那一条边,这条边方向上最远的顶点就是极值点(这里有点像凸包的对踵点算法)。
#include<iostream> #include <cmath> #include <algorithm> #include<vector> using namespace std; #define MAX_N 100000 * 2 typedef long long LL; inline int compare_to_zero(long double x) { if (fabs(x) < 1e-12)return 0; return x < 0 ? -1 : 1; } struct Line { LL a, b, c; } ls[MAX_N], // 输入直线 st[MAX_N]; // 输入st对 long double ans[MAX_N]; bool impossible[MAX_N]; // 优先斜率,截距次之的排序 inline bool cmp(const Line &a, const Line &b) { LL t = a.a*b.b - a.b*b.a; if (t != 0)return t > 0; t = a.c*b.b - a.b*b.c; return t > 0; } struct Point { long double x, y; }; int stack[MAX_N]; Point p[MAX_N]; // 直线交点 inline Point intersection(const Line &la, const Line &lb) { LL a = la.a, b = la.b, c = la.c; LL d = lb.a, e = lb.b, f = lb.c; LL der = a*e - b*d; Point ret; ret.x = c*e - b*f; ret.y = a*f - c*d; ret.x /= der; ret.y /= der; return ret; } // 点是否满足直线的约束 inline bool is_point_fit_line(const Line &l, const Point &p) { return compare_to_zero(l.a * p.x + l.b * p.y - l.c) >= 0; } ///////////////////////////SubMain////////////////////////////////// int main(int argc, char *argv[]) { #ifndef ONLINE_JUDGE freopen("in.txt", "r", stdin); freopen("out.txt", "w", stdout); #endif int n, m; scanf("%d%d", &n, &m); { for (int i = 0; i < n; i++) { scanf("%d%d%d", &ls[i].a, &ls[i].b, &ls[i].c); } sort(ls, ls + n, cmp); int size = 1; for (int i = 1; i < n; i++) { if (ls[i].a* ls[i - 1].b == ls[i - 1].a* ls[i].b) continue; // 排重 ls[size++] = ls[i]; } n = size; LL S, T; if (n == 1) { while (m--) { scanf("%lld%lld", &S, &T); if (S* ls[0].b != T* ls[0].a) puts("IMPOSSIBLE"); // 平行,没有最小值,也就是最小值是无穷小 else printf("%.5f\n", ls[0].c*S*1.0 / ls[0].a); } } else { // 半平面交集 stack[0] = 0; // 一个栈,储存的是凸包边所在直线的下标 p[0].x = 0; // p是最终的交集对应的点,每个点是每条边的最远点,第一个点当然是由x>=0这个半平面和第一个半平面的交点 p[0].y = ls[0].c / (double)ls[0].b; stack[1] = 1; p[1] = intersection(ls[0], ls[1]); // 第二个点 int stack_size = 2; int ls_index; for (ls_index = 2; ls_index < n; ls_index++) { while (stack_size > 1) { int prev = stack[stack_size - 2]; int top = stack[stack_size - 1]; if (is_point_fit_line(ls[top], intersection(ls[ls_index], ls[prev]))) --stack_size; else { break; } } p[stack_size] = intersection(ls[ls_index], ls[stack[stack_size - 1]]); stack[stack_size++] = ls_index; } // 半平面交集求解完毕 for (int i = 0; i < m; i++) { scanf("%d%d", &st[i].a, &st[i].b); st[i].c = i; // 用c表示id impossible[i] = false; } sort(st, st + m, cmp); int st_index; for (st_index = 0; st_index < m; st_index++) { if (st[st_index].a * ls[0].b > st[st_index].b * ls[0].a) impossible[st[st_index].c] = true; // 角度比第一条还小,这条线肯定与半平面交集不相交 else break; } int cur = 0; while (cur < stack_size && st_index < m) { ls_index = stack[stack_size - 1]; if (st[st_index].a * ls[ls_index].b < st[st_index].b * ls[ls_index].a) // 比最后一条还大,还是不相交 { break; } // 角度在半平面交集的范围内,叉积由小增大 while (cur < stack_size) { ls_index = stack[cur]; LL det = st[st_index].a * ls[ls_index].b - st[st_index].b * ls[ls_index].a; // 叉积 if (det >= 0) // 跟这条约束线的方向最相近,在此方向上取最远点作为极值点 { ans[st[st_index].c] = st[st_index].a * p[cur].x + st[st_index].b * p[cur].y; break; } else { cur++; continue; } } st_index++; } while (st_index < m) { impossible[st[st_index].c] = true; st_index++; } for (st_index = 0; st_index < m; st_index++) { if (impossible[st_index]) { puts("IMPOSSIBLE"); } else { printf("%.5f\n", ans[st_index]); } } } } #ifndef ONLINE_JUDGE fclose(stdin); fclose(stdout); system("out.txt"); #endif return 0; } ///////////////////////////End Sub//////////////////////////////////
13915434 | hankcs | 3689 | Accepted | 2888K | 610MS | C++ | 4057B | 2015-02-27 02:51:16 |
其实半平面算法我还找到了一篇文章讲解得不错:http://www.csie.ntnu.edu.tw/~u91029/Half-planeIntersection.html
Reference
http://hi.baidu.com/chenwenwen0210/item/4133c1cbe26308e4994aa050