The MATLAB Blog

Practical Advice for People on the Leading Edge

Linear Programming, the HiGHS Optimization library and MATLAB

What is linear programming?

Linear programming (LP) is a mathematical method used to determine the best possible outcome (such as maximum profit or lowest cost) in a mathematical model whose requirements are represented by linear relationships. It involves optimizing a linear objective function, subject to a set of linear constraints (equalities and inequalities).
Put more formally, we want to find the vector x that minimizes the function
minx{fx}
subject to the constraints
Axb        (inequalityconstraints)Aeqx=beq   (equalityconstraints)lbxub(bound  constraints)
whereAand Aeq are matrices while f, b, beq, lb and ub are vectors.
For me, this makes more sense when looking at a specific example.So, let's consider the classic diet problem as described in [1]. Suppose you have three foods available: Corn, Milk and Bread
Food
Cost / Serving
Vitamin A
Calories
Maximum servings
Corn (C)
$0.18
107
72
10
Milk (M)
$0.23
500
121
10
Bread (B)
$0.05
0
65
10
You also have two nutritional requirements to satisfy
Nutrient
Minimum amount in diet
Maximum amount in diet
Calories
2,000
2,250
Vitamin A
5,000
50,000
and you can have a maximum of 10 portions of any one food type.
How would you choose the foods that minimise the cost while satisfying all of the constraints?
In MATLAB, here's how to solve this using the Problem-based approach, which is the easiest way to get started with solving linear programs in MATLAB. An alternative method would be to use the solver-based approach, which we'll look at later.
% How many servings of each food to eat each day
C = optimvar("corn", LowerBound = 0,UpperBound = 10);
M = optimvar("milk", LowerBound = 0,UpperBound = 10);
B = optimvar("bread",LowerBound = 0,UpperBound = 10);
 
prob = optimproblem;
prob.Objective = C*0.18 + M*0.23 + B*0.05; % The cost to minimise
 
% Add constraints for the Vitamin A requirements
prob.Constraints.c1 = C*107 + M*500 >= 5000;
prob.Constraints.c2 = C*107 + M*500 <= 50000;
 
% Add constraints for the Calories requirements
prob.Constraints.c3 = C*72 + M*121 + B*65 >= 2000;
prob.Constraints.c4 = C*72 + M*121 + B*65 <= 2250;
 
diet = solve(prob)
Solving problem using linprog. Optimal solution found.
diet = struct with fields:
bread: 10 corn: 1.9444 milk: 10
The solve function tells us that it decided to use the linprog function to find the solution, useful information if we ever decide to switch to the-solver based approach in the future.
The result tells us that the cheapest solution is to have 10 portions of bread, 1.944 portions of corn and 10 portions of milk. To see how much this costs, we evaluate the objective function at the solution point.
cost = evaluate(prob.Objective, diet);
fprintf("The optimal diet would cost $%.2f\n",cost)
The optimal diet would cost $3.15
This is a very small problem with only 3 variables; you may even be able to solve it by hand!

Real linear problems are difficult to solve

Real problems tend to be much larger with many thousands of variables and constraints. In addition, there is a class of linear programs called Mixed Integer Linear Programs (MILP) which have the extra constraint that some or all entries of the solution vector be integers.
Solving them can be very difficult! Indeed, there is a collection called MIPLIB 2017 (to which MathWorks is a contributor) that contains problems so hard that even the best solvers in the world can't tackle them. To get your head around this, consider that an easy problem in the MIPLIB 2017 collection is defined as something that can be solved in less than an hour on a machine with 16 cores. That's a lot of compute for something that's defined as 'easy'.
Let's take a look at one such 'easy' problem from MIPLIB 2017 using the latest version of MATLAB, R2024b. The problem is called acc-tight5 and is a basketball scheduling instance that's stored in the MPS file format.
% Download and unzip the MPS file containing the problem
gunzip("https://miplib.zib.de/WebData/instances/acc-tight5.mps.gz")
problem = mpsread("acc-tight5.mps")
Running HiGHS 1.7.0: Copyright (c) 2024 HiGHS under MIT licence terms
problem = struct with fields:
f: [1339×1 double] Aineq: [2775×1339 double] bineq: [2775×1 double] Aeq: [277×1339 double] beq: [277×1 double] lb: [1339×1 double] ub: [1339×1 double] intcon: [1339×1 double] solver: 'intlinprog' options: [1×1 optim.options.Intlinprog]
The mpsread function creates a structure that is suitable for the solver-based approach to solving linear programs and also tells us which solver we should use: intlinprog, since this is a mixed integer problem.
acc_tight5_solution = intlinprog(problem);
Running HiGHS 1.7.0: Copyright (c) 2024 HiGHS under MIT licence terms Coefficient ranges: Matrix [1e+00, 1e+00] Cost [1e+00, 1e+00] Bound [1e+00, 1e+00] RHS [1e+00, 4e+00] Presolving model 3041 rows, 1332 cols, 16057 nonzeros 0s 2454 rows, 1007 cols, 12389 nonzeros 0s 2337 rows, 988 cols, 12345 nonzeros 0s Objective function is integral with scale 1 Solving MIP model with: 2337 rows 988 cols (988 binary, 0 integer, 0 implied int., 0 continuous) 12345 nonzeros Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time 0 0 0 0.00% 0 inf inf 0 0 0 0 0.2s 0 0 0 0.00% 0 inf inf 0 0 3 1635 0.4s 0 0 0 0.00% 0 inf inf 381 10 385 75371 13.6s 0.7% inactive integer columns, restarting Model after restart has 2312 rows, 981 cols (981 bin., 0 int., 0 impl., 0 cont.), and 12154 nonzeros 0 0 0 0.00% 0 inf inf 10 0 0 75371 13.6s 0 0 0 0.00% 0 inf inf 10 6 8 75740 13.7s 0 0 0 0.00% 0 inf inf 384 19 286 135490 25.0s 59 8 22 0.87% 0 inf inf 425 19 1006 175042 30.3s 120 17 50 9.79% 0 inf inf 501 22 2020 212149 35.4s 176 19 74 11.07% 0 inf inf 564 26 2705 251875 40.6s 246 26 104 12.06% 0 inf inf 619 28 3560 289301 45.7s 292 29 122 13.13% 0 inf inf 768 29 4472 327492 50.9s B 311 29 131 13.37% 0 1 100.00% 785 29 4949 339919 52.7s 432 49 173 22.12% 0 1 100.00% 860 9 5728 374411 57.7s 503 52 207 32.03% 0 1 100.00% 890 16 6425 407007 62.7s 578 54 239 35.77% 0 1 100.00% 1361 14 7012 439704 67.8s 682 77 278 36.29% 0 1 100.00% 1279 17 7776 473939 72.9s 763 73 315 39.62% 0 1 100.00% 1597 18 8667 505596 78.0s 858 77 359 40.62% 0 1 100.00% 1702 23 9529 538912 83.1s T 890 0 373 100.00% 0 0 0.00% 1720 27 9794 546054 84.4s Solving report Status Optimal Primal bound 0 Dual bound 0 Gap 0% (tolerance: 0.01%) Solution status feasible 0 (objective) 0 (bound viol.) 2.78645351635e-13 (int. viol.) 0 (row viol.) Timing 84.40 (total) 0.23 (presolve) 0.00 (postsolve) Nodes 891 LP iterations 546171 (total) 49690 (strong br.) 9159 (separation) 129933 (heuristics) Optimal solution found. Intlinprog stopped because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 1e-06. The intcon variables are integer within tolerance, options.ConstraintTolerance = 1e-06.
We can see that intlinprog took 84.4 seconds to solve this problem and that, like the mpsread function, it uses something called HiGHS. We switched to using the HiGHS library by default in MATLAB R2024a for both linprog and intlinprog, while mpsread started using it in R2024b. You can get a clue as to why we switched to HiGHS by asking intlinprog to use its original solver, the one that was used by default until recently.
options = optimoptions("intlinprog",Algorithm="legacy");
Warning: Option "Algorithm" will be removed in a future release. Then intlinprog will use its current default algorithm, "highs".
problem.options = options;
acc_tight5_solution_old = intlinprog(problem);
LP: Optimal objective value is 0.000000. Cut Generation: Applied 81 clique cuts, and 3 zero-half cuts. Lower bound is 0.000000. Cut Generation: Applied 16 clique cuts, and 2 zero-half cuts. Lower bound is 0.000000. Branch and Bound: nodes total num int integer relative explored time (s) solution fval gap (%) 7202 1617.89 1 0.000000e+00 0.000000e+00 Optimal solution found. Intlinprog stopped because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 0. The intcon variables are integer within tolerance, options.IntegerTolerance = 1e-05.
At 1617.9 seconds, the legacy algorithm took rather longer, 20x longer in fact. Now such amazing speed-ups won't always be the case. Sometimes, HiGHS is only a little faster, other times it can even be slower. When considering many problems the transition was worth making.

What is HiGHS?

HiGHS is high performance software for solving large-scale sparse linear programming (LP), mixed-integer programming (MIP) and quadratic programming (QP) models. It is developed in C++11, with interfaces to several other languages.The name HiGHS comes from the initials of its original developers.
Back in June 2024 I had the pleasure of attending the first ever HiGHS workshop at The University of Edinburgh of which MathWorks was one of the sponsors and where MathWorks developer, Franz Wesselmann, was one of the speakers. While there I met Julian Hall and Ivet Galabova who are responsible for three of the initials in HiGHS and are the current maintainers and leaders of the project. MathWorks has been working closely with them for a while now as we worked towards incorporating HiGHS into MATLAB via Optimization Toolbox.
One of the contributions that MathWorks makes back into the HiGHS project is to work directly on the code. MathWorker Franz Wesselmann is now in the top 5 contibutors to the project in terms of numbers of commits.
This is a fruitful collaboration for all concerned and you can expect more HiGHS-related functionality coming to MATLAB in the future. Thanks to everyone in the HiGHS team for welcoming us to their community.

References

|
  • print

コメント

コメントを残すには、ここ をクリックして MathWorks アカウントにサインインするか新しい MathWorks アカウントを作成します。