using JuMP, HiGHS

Let's look at a particular LP \begin{split}\min_{x,y}\quad &-x\\ s.t. \quad &2x + y \leq 1.5\\ & x \geq 0, y \geq 0\end{split}

function mylp()
    model = Model(HiGHS.Optimizer)
    @variable(model, 0 <= x)
    @variable(model, 0 <= y)
    @objective(model, Min, -x)
    @constraint(model, 2x + y <= 1.5)
    optimize!(model)
    return (
        status = termination_status(model),
        objval = objective_value(model),
        sol = (value.(x), value.(y))
    )
end 

mylp()

Running HiGHS 1.9.0 (git hash: 66f735e60): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 2e+00]
  Cost   [1e+00, 1e+00]
  Bound  [0e+00, 0e+00]
  RHS    [2e+00, 2e+00]
Presolving model
1 rows, 1 cols, 1 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve : Reductions: rows 0(-1); columns 0(-2); elements 0(-2) - Reduced to empty
Solving the original LP from the solution after postsolve
Model status        : Optimal
Objective value     : -7.5000000000e-01
Relative P-D gap    :  0.0000000000e+00
HiGHS run time      :          0.00

(status = MathOptInterface.OPTIMAL, objval = -0.75, sol = (0.75, 0.0))

Form 1

linprog(c, A, lb, ub, l, u)

\begin{split}\min_{x}\quad &c^Tx\\ s.t.\quad &lb \leq Ax \leq ub\\ &l \leq x \leq u\\\end{split}

using JuMP

_fix_length(x,N) = length(x) == 1 ? x*ones(N) : x 
function linprog(c, A, lb::Union{Real,Vector}, ub, l, u, solver)
    N = length(c)
    model = Model(solver)
    l = _fix_length(l,N)
    u = _fix_length(u,N)
    @variable(model, l[i] .<= x[i=1:N] .<= u[i])
    @objective(model, Min, c' * x)
    @constraint(model, lb .<= A*x .<= ub)
    optimize!(model)
    return (
        status = termination_status(model),
        objval = objective_value(model),
        sol = value.(x)
    )
end

linprog (generic function with 1 method)

Let's look at a particular LP \begin{split}\min_{x,y}\quad &-x\\ s.t. \quad &2x + y \leq 1.5\\ & x \geq 0, y \geq 0\end{split}

sol = linprog([-1,0],[2 1],-Inf,1.5,0,Inf,HiGHS.Optimizer)

Running HiGHS 1.9.0 (git hash: 66f735e60): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 2e+00]
  Cost   [1e+00, 1e+00]
  Bound  [0e+00, 0e+00]
  RHS    [2e+00, 2e+00]
Presolving model
1 rows, 1 cols, 1 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve : Reductions: rows 0(-1); columns 0(-2); elements 0(-2) - Reduced to empty
Solving the original LP from the solution after postsolve
Model status        : Optimal
Objective value     : -7.5000000000e-01
Relative P-D gap    :  0.0000000000e+00
HiGHS run time      :          0.00

(status = MathOptInterface.OPTIMAL, objval = -0.75, sol = [0.75, 0.0])

sol.objval

-0.75

sol.sol

2-element Vector{Float64}:
 0.75
 0.0

sol.status

OPTIMAL::TerminationStatusCode = 1

Form 2

linprog(c, A, sense, b)

\begin{split}\min_{x}\quad &c^Tx\\ s.t. \quad &a_i^Tx \text{ sense}_i \, b_i \forall\,\, i\\ &0 \leq x\\\end{split}

""" This function has an implicit x >= 0 constraint """
function linprog(c, A, sense, b, solver)
    N = length(c)
    model = Model(solver)
    set_silent(model)
    @variable(model, x[i=1:N] >= 0)
    @objective(model, Min, c' * x)
    eq_rows, ge_rows, le_rows = sense .== '=', sense .== '>', sense .== '<'
    @constraint(model, A[eq_rows, :] * x .== b[eq_rows])
    @constraint(model, A[ge_rows, :] * x .>= b[ge_rows])
    @constraint(model, A[le_rows, :] * x .<= b[le_rows])
    optimize!(model)
    return (
        status = termination_status(model),
        objval = objective_value(model),
        sol = value.(x)
    )
end

linprog

sol = linprog([-1,0],[2 1],['<'],[1.5],HiGHS.Optimizer)

(status = MathOptInterface.OPTIMAL, objval = -0.75, sol = [0.75, 0.0])

sol = linprog([-1,0],[2 1],['='],[1.5],HiGHS.Optimizer)

(status = MathOptInterface.OPTIMAL, objval = -0.75, sol = [0.75, 0.0])

sol = linprog([-1,0],[2 1],['>'],[1.5],HiGHS.Optimizer)

(status = MathOptInterface.DUAL_INFEASIBLE, objval = -0.5, sol = [0.5, 0.0])

sol = linprog([-1,0],[2 1; 1 0],['<','<'],[1.5; -0.25],HiGHS.Optimizer)

(status = MathOptInterface.INFEASIBLE, objval = 0.0, sol = [5.0e-324, 1.7e-322])

sol.objval