{"id":2834,"date":"2024-11-04T11:11:23","date_gmt":"2024-11-04T16:11:23","guid":{"rendered":"https:\/\/blogs.mathworks.com\/matlab\/?p=2834"},"modified":"2024-11-04T11:11:23","modified_gmt":"2024-11-04T16:11:23","slug":"linear-programming-the-highs-optimization-library-and-matlab","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/matlab\/2024\/11\/04\/linear-programming-the-highs-optimization-library-and-matlab\/","title":{"rendered":"Linear Programming, the HiGHS Optimization library and MATLAB"},"content":{"rendered":"

What is linear <\/span>programming<\/span>?<\/span><\/h4>
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). <\/span><\/div>
Put more formally, we want to find the vector x that minimizes the function<\/span><\/div>
min<\/mi><\/mrow>x<\/mi><\/mrow><\/munder>{<\/mo>f<\/mi><\/mrow>\u22a4<\/mo><\/mrow><\/msup>x<\/mi><\/mrow>}<\/mo><\/mrow><\/mrow><\/math><\/span><\/div>
subject to the constraints <\/span><\/div>
Ax<\/mi>\u2264<\/mo>b<\/mi>\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009<\/mtext>(<\/mo>inequality<\/mi>\u2009<\/mtext>constraints<\/mi><\/mrow>)<\/mo><\/mrow><\/mrow><\/mtd><\/mtr>A<\/mi><\/mrow>eq<\/mi><\/mrow><\/msub>x<\/mi>=<\/mo>b<\/mi><\/mrow>eq<\/mi><\/mrow><\/msub>\u2009\u2009\u2009<\/mtext>(<\/mo>equality<\/mi>\u2009<\/mtext>constraints<\/mi><\/mrow>)<\/mo><\/mrow><\/mrow><\/mtd><\/mtr>lb<\/mi>\u2264<\/mo>x<\/mi>\u2264<\/mo>ub<\/mi>\u2009<\/mtext>(<\/mo>bound<\/mi>\u2009\u2009<\/mtext>constraints<\/mi><\/mrow>)<\/mo><\/mrow><\/mrow><\/mtd><\/mtr><\/mtable><\/mrow><\/math><\/span><\/div>
where<\/span>A<\/span>and <\/span>A<\/mi><\/mrow>eq<\/mi><\/mrow><\/msub><\/mrow><\/math><\/span> are matrices while <\/span>f<\/span>, <\/span>b<\/span>, <\/span>b<\/mi><\/mrow>eq<\/mi><\/mrow><\/msub><\/mrow><\/math><\/span>, <\/span>lb<\/span> and <\/span>ub<\/span> are vectors.<\/span><\/div>
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<\/span><\/div>
Food<\/span><\/div><\/td>
Cost \/ Serving<\/span><\/div><\/td>
Vitamin A<\/span><\/div><\/td>
Calories<\/span><\/div><\/td>
Maximum servings<\/span><\/div><\/td><\/tr>
Corn (C)<\/span><\/div><\/td>
$0.18<\/span><\/div><\/td>
107<\/span><\/div><\/td>
72<\/span><\/div><\/td>
10<\/span><\/div><\/td><\/tr>
Milk (M)<\/span><\/div><\/td>
$0.23<\/span><\/div><\/td>
500<\/span><\/div><\/td>
121<\/span><\/div><\/td>
10<\/span><\/div><\/td><\/tr>
Bread (B)<\/span><\/div><\/td>
$0.05<\/span><\/div><\/td>
0<\/span><\/div><\/td>
65<\/span><\/div><\/td>
10<\/span><\/div><\/td><\/tr><\/table>
You also have two nutritional requirements to satisfy<\/span><\/div>
Nutrient<\/span><\/div><\/td>
Minimum amount in diet<\/span><\/div><\/td>
Maximum amount in diet<\/span><\/div><\/td><\/tr>
Calories<\/span><\/div><\/td>
2,000<\/span><\/div><\/td>
2,250<\/span><\/div><\/td><\/tr>
Vitamin A<\/span><\/div><\/td>
5,000<\/span><\/div><\/td>
50,000<\/span><\/div><\/td><\/tr><\/table>
and you can have a maximum of 10 portions of any one food type. <\/span><\/div>
How would you choose the foods that minimise the cost while satisfying all of the constraints? <\/span><\/div>
In MATLAB, here's how to solve this using the <\/span>Problem-based approach<\/span><\/a>, which is the easiest way to get started with solving linear programs in MATLAB. An alternative method would be to use the <\/span>solver-based approach<\/span><\/a>, which we'll look at later.<\/span><\/div>
% How many servings of each food to eat each day<\/span><\/span><\/div><\/div>
C = optimvar(<\/span>\"corn\"<\/span>, LowerBound = 0,UpperBound = 10);<\/span><\/span><\/div><\/div>
M = optimvar(<\/span>\"milk\"<\/span>, LowerBound = 0,UpperBound = 10);<\/span><\/span><\/div><\/div>
B = optimvar(<\/span>\"bread\"<\/span>,LowerBound = 0,UpperBound = 10);<\/span><\/span><\/div><\/div>
 <\/div><\/div>
prob = optimproblem;<\/span><\/span><\/div><\/div>
prob.Objective = C*0.18 + M*0.23 + B*0.05; <\/span>% The cost to minimise<\/span><\/span><\/div><\/div>
 <\/div><\/div>
% Add constraints for the Vitamin A requirements<\/span><\/span><\/div><\/div>
prob.Constraints.c1 = C*107 + M*500 >= 5000;<\/span><\/span><\/div><\/div>
prob.Constraints.c2 = C*107 + M*500 <= 50000;<\/span><\/span><\/div><\/div>
 <\/div><\/div>
% Add constraints for the Calories requirements<\/span><\/span><\/div><\/div>
prob.Constraints.c3 = C*72 + M*121 + B*65 >= 2000;<\/span><\/span><\/div><\/div>
prob.Constraints.c4 = C*72 + M*121 + B*65 <= 2250;<\/span><\/span><\/div><\/div>
 <\/div><\/div>
diet = solve(prob)<\/span><\/span><\/div>
Solving problem using linprog.\n\nOptimal solution found.<\/div><\/div>
diet = struct with fields:<\/span><\/span><\/div>
bread: 10\n corn: 1.9444\n milk: 10\n<\/div><\/div><\/div><\/div><\/div><\/div>
The <\/span>solve<\/span><\/a> function tells us that it decided to use the <\/span>linprog<\/span><\/a> function to find the solution, useful information if we ever decide to switch to the-solver based approach in the future. <\/span><\/div>
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.<\/span><\/div>
cost = evaluate(prob.Objective, diet);<\/span><\/span><\/div><\/div>
fprintf(<\/span>\"The optimal diet would cost $%.2f\\n\"<\/span>,cost)<\/span><\/span><\/div>
The optimal diet would cost $3.15<\/div><\/div><\/div><\/div><\/div>
This is a very small problem with only 3 variables; you may even be able to solve it by hand! <\/span><\/div>

Real linear problems are difficult to solve<\/span><\/h4>
Real problems tend to be much larger with many thousands of variables and constraints. In addition, there is a class of linear programs called <\/span>Mixed Integer Linear Programs (MILP) <\/span><\/a>which have the extra constraint that some or all entries of the solution vector be integers. <\/span><\/div>
Solving them can be very difficult! Indeed, there is a <\/span>collection called MIPLIB 2017<\/span><\/a> (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 <\/span>easy<\/span> 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'.<\/span><\/div>
Let's take a look at one such 'easy' problem from MIPLIB 2017 using the latest version of MATLAB, R2024b. The <\/span>problem is called acc-tight5<\/span><\/a> and is a basketball scheduling instance that's stored in the <\/span>MPS file format<\/span><\/a>.<\/span><\/div>
% Download and unzip the MPS file containing the problem<\/span><\/span><\/div><\/div>
gunzip(<\/span>\"https:\/\/miplib.zib.de\/WebData\/instances\/acc-tight5.mps.gz\"<\/span>)<\/span><\/span><\/div><\/div>
problem = mpsread(<\/span>\"acc-tight5.mps\"<\/span>)<\/span><\/span><\/div>
Running HiGHS 1.7.0: Copyright (c) 2024 HiGHS under MIT licence terms<\/div><\/div>
problem = struct with fields:<\/span><\/span><\/div>
f: [1339\u00d71 double]\n Aineq: [2775\u00d71339 double]\n bineq: [2775\u00d71 double]\n Aeq: [277\u00d71339 double]\n beq: [277\u00d71 double]\n lb: [1339\u00d71 double]\n ub: [1339\u00d71 double]\n intcon: [1339\u00d71 double]\n solver: 'intlinprog'\n options: [1\u00d71 optim.options.Intlinprog]\n<\/div><\/div><\/div><\/div><\/div><\/div>
The <\/span>mpsread<\/span><\/a> 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: <\/span>intlinprog<\/span><\/a>, since this is a mixed integer problem. <\/span><\/div>
acc_tight5_solution = intlinprog(problem);<\/span><\/span><\/div>
Running HiGHS 1.7.0: Copyright (c) 2024 HiGHS under MIT licence terms\nCoefficient ranges:\n Matrix [1e+00, 1e+00]\n Cost [1e+00, 1e+00]\n Bound [1e+00, 1e+00]\n RHS [1e+00, 4e+00]\nPresolving model\n3041 rows, 1332 cols, 16057 nonzeros 0s\n2454 rows, 1007 cols, 12389 nonzeros 0s\n2337 rows, 988 cols, 12345 nonzeros 0s\nObjective function is integral with scale 1\n\nSolving MIP model with:\n 2337 rows\n 988 cols (988 binary, 0 integer, 0 implied int., 0 continuous)\n 12345 nonzeros\n\n Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work \n Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time\n\n 0 0 0 0.00% 0 inf inf 0 0 0 0 0.2s\n 0 0 0 0.00% 0 inf inf 0 0 3 1635 0.4s\n 0 0 0 0.00% 0 inf inf 381 10 385 75371 13.6s\n\n0.7% inactive integer columns, restarting\nModel after restart has 2312 rows, 981 cols (981 bin., 0 int., 0 impl., 0 cont.), and 12154 nonzeros\n\n 0 0 0 0.00% 0 inf inf 10 0 0 75371 13.6s\n 0 0 0 0.00% 0 inf inf 10 6 8 75740 13.7s\n 0 0 0 0.00% 0 inf inf 384 19 286 135490 25.0s\n 59 8 22 0.87% 0 inf inf 425 19 1006 175042 30.3s\n 120 17 50 9.79% 0 inf inf 501 22 2020 212149 35.4s\n 176 19 74 11.07% 0 inf inf 564 26 2705 251875 40.6s\n 246 26 104 12.06% 0 inf inf 619 28 3560 289301 45.7s\n 292 29 122 13.13% 0 inf inf 768 29 4472 327492 50.9s\n B 311 29 131 13.37% 0 1 100.00% 785 29 4949 339919 52.7s\n 432 49 173 22.12% 0 1 100.00% 860 9 5728 374411 57.7s\n 503 52 207 32.03% 0 1 100.00% 890 16 6425 407007 62.7s\n 578 54 239 35.77% 0 1 100.00% 1361 14 7012 439704 67.8s\n 682 77 278 36.29% 0 1 100.00% 1279 17 7776 473939 72.9s\n 763 73 315 39.62% 0 1 100.00% 1597 18 8667 505596 78.0s\n 858 77 359 40.62% 0 1 100.00% 1702 23 9529 538912 83.1s\n T 890 0 373 100.00% 0 0 0.00% 1720 27 9794 546054 84.4s\n\nSolving report\n Status Optimal\n Primal bound 0\n Dual bound 0\n Gap 0% (tolerance: 0.01%)\n Solution status feasible\n 0 (objective)\n 0 (bound viol.)\n 2.78645351635e-13 (int. viol.)\n 0 (row viol.)\n Timing 84.40 (total)\n 0.23 (presolve)\n 0.00 (postsolve)\n Nodes 891\n LP iterations 546171 (total)\n 49690 (strong br.)\n 9159 (separation)\n 129933 (heuristics)\n\nOptimal solution found.\n\nIntlinprog 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.<\/div><\/div><\/div><\/div><\/div>
We can see that <\/span>intlinprog<\/span> took 84.4 seconds to solve this problem and that, like the <\/span>mpsread<\/span> function, it uses something called HiGHS. We switched to using the HiGHS library by default in MATLAB R2024a for both <\/span>linprog<\/span> and <\/span>intlinprog, <\/span>while <\/span>mpsread<\/span> started using it in R2024b. You can get a clue as to why we switched to HiGHS by asking <\/span>intlinprog<\/span> to use its original solver, the one that was used by default until recently.<\/span><\/div>
options = optimoptions(<\/span>\"intlinprog\"<\/span>,Algorithm=<\/span>\"legacy\"<\/span>);<\/span><\/span><\/div>
Warning: Option \"Algorithm\" will be removed in a future release. Then intlinprog will use its current default algorithm, \"highs\".<\/div>
<\/div><\/div><\/div><\/div><\/div>
problem.options = options;<\/span><\/span><\/div><\/div>
acc_tight5_solution_old = intlinprog(problem);<\/span><\/span><\/div>
LP: Optimal objective value is 0.000000. \n\nCut Generation: Applied 81 clique cuts, \n and 3 zero-half cuts. \n Lower bound is 0.000000. \n\nCut Generation: Applied 16 clique cuts, \n and 2 zero-half cuts. \n Lower bound is 0.000000. \n\nBranch and Bound:\n\n nodes total num int integer relative \nexplored time (s) solution fval gap (%) \n 7202 1617.89 1 0.000000e+00 0.000000e+00 \n\nOptimal solution found.\n\nIntlinprog 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.<\/div><\/div><\/div><\/div><\/div>
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. <\/span><\/div>

What is HiGHS?<\/span><\/h4>
HiGHS<\/span><\/a> 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.<\/span><\/div>
Back in June 2024 I had the pleasure of attending the first ever<\/span> HiGHS workshop<\/span><\/a> 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 <\/span>Julian Hall<\/span><\/a> and <\/span>Ivet Galabova<\/span><\/a> 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 <\/span>incorporating HiGHS into MATLAB<\/span><\/a> via Optimization Toolbox. <\/span><\/div>
\"\"<\/img><\/div>
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<\/span> top 5 contibutors to the project<\/span><\/a> in terms of numbers of commits. <\/span><\/div>
\"\"<\/img><\/a><\/div>
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. <\/span><\/div>

References<\/span><\/h2>
[1] <\/span>The Diet Problem: A WWW-based Interactive Case Study in Linear Programming<\/span><\/a>. Joseph Czyzyk and Timothy J. Wisniewski <\/span><\/div>
<\/div>\n<\/div>