Linear Programming
Mathwrist Code Example Series

Copyright ©Mathwrist LLC 2023
(January 20, 2023)
Abstract

This document gives some code examples on how to use the linear programming (LP) features in Mathwrist’s C++ Numerical Programming Library (NPL).

1 Simplex Method

This example illustrates solving the following LP problem using simplex method.

arg⁡minx,y-4⁢x-y, s.t.
x+y≤50
3⁢x+y≤90
x≥0,y≥0
1   // 2 general linear constraints, 2 unknowns
2   LinearProg::Constraints constr(2, 2);
4   // x + y <= 50
5   constr.A(0, 0) = 1;
6   constr.A(0, 1) = 1;
7   constr.general_bounds[0] = Bound::upper(50);
9   // 3x + y <= 90
10   constr.A(1, 0) = 3;
11   constr.A(1, 1) = 1;
12   constr.general_bounds[1] = Bound::upper(90);
14   // simple bounds, x >= 0, y >= 0
15   constr.simple_bounds[0] = Bound::lower(0);
16   constr.simple_bounds[1] = Bound::lower(0);
18   // cost -4x - y
19   Matrix cost(2, 1);
20   cost(0, 0) = -4;
21   cost(1, 0) = -1;
23   // Initial guess x = 0;
24   Matrix x(2, 1);
26   Simplex simplex(constr);
27   double z = simplex(cost, x);
29   assert(std::abs(x(0, 0) - 30) < 1.e-10);
30   assert(std::abs(x(1, 0)) < 1.e-10);
31   assert(std::abs(z + 120) < 1.e-10);

2 LP Active Set Method

This example illustrates solving the following LP problem using LP active-set method.

arg⁡minx,y⁡200⁢x+500⁢y, s.t.
x+2⁢y≥10
3⁢x+4⁢y≤24
x≥0,y≥0
1 // Linear constraints:
2 // x  + 2y >= 10
3 // 3x + 4y <= 24
5 // Matrix view from constraint coefficients
6 const double _A[] = { 1, 2, 3, 4 };
7 Matrix::ViewConst A(2, 2, _A);
9 // constraint bounds;
10 ftl::Buffer<Bound> b(2);
11 b[0] = Bound::lower(10);
12 b[1] = Bound::upper(24);
14 // Create constraint specification from coefficient matrix
15 // and bounds.
16 LinearProg::Constraints constr(A, b);
18 // Simple bounds, x >= 0, y >= 0
19 constr.simple_bounds[0] = Bound::lower(0);
20 constr.simple_bounds[1] = Bound::lower(0);
22 // cost : 200 x + 500 y
23 Matrix c(2, 1);
24 c(0, 0) = 200;
25 c(1, 0) = 500;
27 LP_ActiveSet lp(constr);
29 // Initial guess x = 0;
30 Matrix x(2, 1);
32 lp.set_crash_start_radius(10);
33 double z = lp(c, x);
35 assert(std::abs(x(0, 0) - 4) < 1.e-10);
36 assert(std::abs(x(1, 0) - 3) < 1.e-10);
37 assert(std::abs(z - 2300) < 1.e-10);