使用lpsolve解决线性规划问题
Solve LP problem in lpsolve
一、引言 Introduction
通过一个简单例子来介绍lpsolve求解线性规划问题的方法。
假若农民有75亩地，他打算种上两种农作物：小麦和大麦。为了种植这些农作物，农民在种子和化肥等的开销分别为：小麦每亩需要$120，大麦每亩为$210。这个农民可支出的钱有$15000。但是当这些农作物丰收后，需要有仓库来存储，以便在行情最好的时候将这些农作物买出。农民的仓库可存储4000蒲式耳(计量谷物等的容量单位，在英国等于36.368升，在美国等于35.238升)。每亩小麦平均产量为110蒲式耳，每亩大麦平均产量为30蒲式耳。若每蒲式耳小麦的净利润为$1.30，每蒲式耳大麦的净利润为$2.00，农民怎么样来种植这75亩地能获得最大利益？
二、建立数学模型 Formulation of an lp problem
首先，我们求出目标函数（Objective），即利润；然后找出约束条件，并画出图形；最后，我们通过图形和一些简单计算得到最优解。
设小麦的种植面积为x，大麦的种植面积为y。
则P=(110)(1.30)x + (30)(2.00)y = 143x + 60y 为目标函数。当其取最大值时也表示农民将获得最大利润。还有如下三个约束不等式，约束分别来自可用支出的限制、存储容量限制和种植面积限制。分别表示如下：
120x + 210y <= 15000
110x + 20y <= 4000
x + y <= 75
严格来说，还有两个不等式的约束，即非负约束来自实际情况，因为农民不可能在负数的地上种植庄稼。即：x >= 0, y >= 0。
得数学模型如下所示：
max(143x + 60y)
s.t.
120x + 210y <= 15000
110x + 30y <= 4000
x + y <= 75
x >= 0
y >= 0
三、图解法
非负约束表明我们只需要考虑X-Y平面中的第一象限即可。x + y <= 75表示以直线x+y=75为界的左下方区域。作图如下所示：
加上另外两个约束条件后，如下图所示：
黑色区域为可行域，即在可行域中任意一点都满足约束条件。
沿目标函数梯度方向作等值线，如下图所示：
等值线与黑色区域相交均为可行解。等值线越向右移动，则目标函数值越大。最优解就是取得最大值且与黑色区域仍有交线的那条等值线。
通过图示，可以清楚地看出目标函数P的最大值在可行域的多边形顶点上取得，坐标值大约为（22, 53）。即直线x+y=75与直线110*x+30*y=4000的交点。这个点是可行域多边形的角点。这不是偶然现象，lpsolve中使用的单纯形算法求得的最优解也是在角点处取得。
定理：对线性规划问题，若可行域有界且存在最优解，则目标函数必可在其可行域的某个顶点达到最优值。
四、使用C/C++建立模型
上述模型在C中建立如下：
1: //------------------------------------------------------------------------------
2: // Copyright (c) 2012 eryar All Rights Reserved.
3: //
4: // File : Main.cpp
5: // Author : eryar@163.com
6: // Date : 2012-9-9 17:11
7: // Version : 0.1v
8: //
9: // Description : lpsolve test program.
10: //
11: //==============================================================================
12:
13: #include <iostream>
14: using namespace std;
15:
16: #include "lp_lib.h"
17:
18: #pragma comment(lib, "liblpsolve55.lib")
19:
20: int demo(void);
21:
22: int main(int argc, char* argv[])
23: {
24: return demo();
25: }
26:
27: int demo( void )
28: {
29: lprec* lp;
30: int Ncol = 0;
31: int *colno = NULL;
32: int j = 0;
33: int ret = 0;
34: REAL* row = NULL;
35:
36: // We will build the model row by row,
37: // So we start with creating a model with 0 rows and 2 columns.
38:
39: // There are two bariables in the model.
40: Ncol = 2;
41:
42: lp = make_lp(0, Ncol);
43:
44: if (lp == NULL)
45: {
46: cout<<"Unable to create new LP model!\n"<<endl;
47:
48: ret = 1;
49: }
50:
51: if (ret == 0)
52: {
53: // Let us name our bariables. Not required, but can be useful for debugging.
54: set_col_name(lp, 1, "x");
55: set_col_name(lp, 2, "y");
56:
57: // Create space large enough for one row.
58: colno = (int *)malloc(Ncol * sizeof(*colno));
59: row = (REAL*)malloc(Ncol * sizeof(*row));
60:
61: // malloc memory failed.
62: if ((colno == NULL) || row == NULL)
63: {
64: ret = 2;
65: }
66: }
67:
68: if (ret == 0)
69: {
70: // Makes building the model faster if it is done rows by row.
71: set_add_rowmode(lp, TRUE);
72:
73: // Construct first row (120 x + 210 y <= 15000).
74: j = 0;
75:
76: // First column.
77: colno[j] = 1;
78: row[j++] = 120;
79:
80: // Second column.
81: colno[j] = 2;
82: row[j++] = 210;
83:
84: // Add the row to lpsolve.
85: if (!add_constraintex(lp, j, row, colno, LE, 15000))
86: {
87: ret = 3;
88: }
89: }
90:
91: // Construct second row (110x + 30y <= 4000).
92: if (ret == 0)
93: {
94: j = 0;
95:
96: // First column.
97: colno[j] = 1;
98: row[j++] = 110;
99:
100: // Second column.
101: colno[j] = 2;
102: row[j++] = 30;
103:
104: // Add the row to lpsolve.
105: if (!add_constraintex(lp, j, row, colno, LE, 4000))
106: {
107: ret = 3;
108: }
109: }
110:
111: // Construct third row (x + y <= 75).
112: if (ret == 0)
113: {
114: j = 0;
115:
116: // First column.
117: colno[j] = 1;
118: row[j++] = 1;
119:
120: // Second column.
121: colno[j] = 2;
122: row[j++] = 1;
123:
124: // Add the row the lpsolve.
125: if (!add_constraintex(lp, j, row, colno, LE, 75))
126: {
127: ret = 3;
128: }
129: }
130:
131: // Set the objective function (143x + 60y).
132: if (ret == 0)
133: {
134: // Row mode should be truned off again when done building the model.
135: set_add_rowmode(lp, FALSE);
136:
137: // Set the objective function (143x + 60y).
138: j = 0;
139:
140: // First column.
141: colno[j] = 1;
142: row[j++] = 143;
143:
144: // Second column.
145: colno[j] = 2;
146: row[j++] = 60;
147:
148: // Set the objective in lpsolve.
149: if (!set_obj_fnex(lp, j, row, colno))
150: {
151: ret = 4;
152: }
153: }
154:
155: if (ret == 0)
156: {
157: // Set the object direction to maximize.
158: set_maxim(lp);
159:
160: // Just out of curioucity, now show the model in lp format on screen.
161: // This only works if this is a console application. If not, use
162: // wirte_lp and a file name.
163: write_LP(lp, stdout);
164:
165: // I only want to see important messages on screen while solving.
166: set_verbose(lp, IMPORTANT);
167:
168: // Now let lpsolve calculate a solution.
169: ret = solve(lp);
170:
171: if (ret == OPTIMAL)
172: {
173: ret = 0;
174: }
175: else
176: {
177: ret = 5;
178: }
179: }
180:
181: // A solution is calculated, now lets get some results.
182: if (ret == 0)
183: {
184: // Objective value.
185: cout<<"Objective value: "<<get_objective(lp)<<endl;
186:
187: // Variable values.
188: get_variables(lp, row);
189:
190: for (j = 0; j < Ncol; j++)
191: {
192: cout<<get_col_name(lp, j+1)<<"="<<row[j]<<endl;
193: }
194:
195: // We are done now.
196: }
197:
198: // Free allocated memory.
199: if (row != NULL)
200: {
201: free(row);
202: }
203:
204: if (colno != NULL)
205: {
206: free(colno);
207: }
208:
209: // Clean up such that all used memory by lpsolve is freed.
210: if (lp != NULL)
211: {
212: delete_lp(lp);
213: }
214:
215: return ret;
216: }
217:
计算结果如下所示：
五、结论
通过这个简单例子，对lpsolve的使用有了个初步认识。若想进一步，可以参考其文档。
PDF Version: Solve LP problem in lpsolve