Solving a linear system of equalities and inequalities

Specific Domains Optimization (Mathematical)
1

I have a mixed system of equalities and inequalities, for example:

x + 2y +3z = 10
2x + 4y + 10z = 20
4x + y + z < 10

where unknown variables x, y, z are all real-valued.

How can get a solution (or the range of all feasible solutions) for this system in Julia? If the system were all equalities, I can format it as a matrix multiplication and use the \ operator for the solution. I understand this can be converted into a linear constrained optimisation problem, but I do not know how. What packages should I use?

Please help

2

with JuMP.jl
although you will get one solution only, not a range

1 Like
3

Thanks! I’ve had a look at JuMP.jl, but to formulate this as an optimisation problem I need to define an objective function to maximise/minimise (and the constraints themselves). Is there a resource that can show me how to do this for my system?

4

you don’t have to give an objective, you will just get a (random ?) feasible solution

3 Likes
5

here is your simple example with JuMP

using JuMP, Clp

m = Model(Clp.Optimizer)

@variables m begin
     x
     y
     z
end

@constraints m begin
      x + 2y +  3z == 10
     2x + 4y + 10z == 20
     4x +  y +   z <= 10
 end

optimize!(m)

value.((x,y,z)) # (0.0, 5.0, 0.0)
9 Likes
6

Thanks! I will give this a try. Is there any way to return the range of all feasible solutions? e.g. in this example the inequality constraint results in y > 30/7. This is obviously simple in this example, but my system has a very large number of equalities and inequalities, and it would be useful to inspect the properties of the space of feasible solutions

7

I do not think this is possible even in general. The range of valid values for a variable y depends on the value of the variable x and z. Unless you mean the first y value for which exists a valid solution and the last y value for which a valid solution (for whichever other variables values). There is something called “sensibility analysis” that goes in the direction you want, but I never dabbled on it, and I do not know its limitations or if it is implemented currently in JuMP.

8

you can get the matrix representation of your problem with this function if it helps
(same problem as before)

julia> Asparse, lb, ub, _ = JuMP._std_matrix(m);

julia> collect(Asparse)
3×6 Array{Float64,2}:
 1.0  2.0   3.0  -1.0   0.0   0.0
 2.0  4.0  10.0   0.0  -1.0   0.0
 4.0  1.0   1.0   0.0   0.0  -1.0

julia> lb
6-element Array{Float64,1}:
 -Inf
 -Inf
 -Inf
  10.0
  20.0
 -Inf

julia> ub
6-element Array{Float64,1}:
 Inf
 Inf
 Inf
 10.0
 20.0
 10.0
9

Oh that’s very interesting, thanks for the help!

10

All feasible solutions are described by a polytope (or the relative interior in case of strict inequalities). So a software package which can compute polytopes, e.g., Polymake from a list of equalities and inequalities can solve your problem.

5 Likes
11

If I try this example with a large system with 100+ variables and as many constraints, I get the error “Primal solution not available”. But I know a solution exists. Does the Clp optimizer not behave well with very large dimensions?

12

Perhaps see if youf problem fits into an Linear Complementarity Problem (LCP) or MCP . If so, https://github.com/odow/PATH.jl is solid. You can use it with Jump or without, and there is an academic license if required for large problems.

3 Likes
13

As others have mentioned, the problem you posted can be formulated as a feasibility LP. Here is how to use the LazySets.jl API for such problem.

The set p constructed below represents the entire solution set of the problem.

using LazySets, ModelingToolkit, Plots

# define modeling variables
vars = @variables x y z

# input the polyhedron in constraint representation
p = HPolyhedron([x + 2y + 3z == 10,
                 2x + 4y + 10z == 20,
                 4x + y + z < 10], vars)

# remove redundant inequalities
remove_redundant_constraints!(p)

# print the list of constraints
constraints_list(p)
5-element Array{HalfSpace{Float64,Array{Float64,1}},1}:
 HalfSpace{Float64,Array{Float64,1}}([1.0, 2.0, 3.0], 10.0)
 HalfSpace{Float64,Array{Float64,1}}([-1.0, -2.0, -3.0], -10.0)
 HalfSpace{Float64,Array{Float64,1}}([2.0, 4.0, 10.0], 20.0)
 HalfSpace{Float64,Array{Float64,1}}([-2.0, -4.0, -10.0], -20.0)
 HalfSpace{Float64,Array{Float64,1}}([4.0, 1.0, 1.0], 10.0)

# get the matrix representation Ax <= b  (where <= is taken component-wise)
tosimplehrep(p)
([1.0 2.0 3.0; -1.0 -2.0 -3.0; … ; -2.0 -4.0 -10.0; 4.0 1.0 1.0], [10.0, -10.0, 20.0, -20.0, 10.0])

There are many methods applicable to polyhedra in constraint representation. For instance,

# check emptiness
isempty(p)
false

# query an element
z = an_element(p)
3-element Array{Float64,1}:
 1.4285714285714286
 4.285714285714286
 0.0

# check membership
z ∈ p
true

-z ∈ p
false

# solution found in another comment using Clp optimizer
[0.0, 5.0, 0.0] ∈ p
true

# check inclusion
Ball2(z, 0.1) ⊆ p
false

# lazy projection
q = Projection(p, 1:2)
LinearMap{Float64,HPolyhedron{Float64,Array{Float64,1}},Float64,SparseArrays.SparseMatrixCSC{Float64,Int64}}(
  [1, 1]  =  1.0
  [2, 2]  =  1.0, HPolyhedron{Float64,Array{Float64,1}}(HalfSpace{Float64,Array{Float64,1}}[HalfSpace{Float64,Array{Float64,1}}([1.0, 2.0, 3.0], 10.0), HalfSpace{Float64,Array{Float64,1}}([-1.0, -2.0, -3.0], -10.0), HalfSpace{Float64,Array{Float64,1}}([2.0, 4.0, 10.0], 20.0), HalfSpace{Float64,Array{Float64,1}}([-2.0, -4.0, -10.0], -20.0), HalfSpace{Float64,Array{Float64,1}}([4.0, 1.0, 1.0], 10.0)]))

I stop there, but please note that to work in high dimensions in an efficient manner it is often a good idea to make lazy operations as above (in the call to Projection). The code that I showed is not limited to 3D, you can try with @variables x[1:100].

3 Likes
14

I get the error “Primal solution not available”. But I know a solution exists since I have a different algorithm that just lowers the RHS of the inequalities (so they become equalities, but this is obviously not rigorous) by some factor and solve the system as a multidimensional root. Does the Clp optimizer not behave well with very large dimensions?

You can check the status with termination_status(model) and primal_status(model). See https://jump.dev/JuMP.jl/stable/solutions/#Obtaining-solutions-1

2 Likes
15

That’s incredible. Thank you for showing me this example!
I have a question that is less implementation now and more maths, since I am lacking in my undergraduate studies in this field. In the Clp.Optimizer example above, the inequality has been implemented as a constraint with the <= operator rather than a strict inequality <. How does this affect the calculations, would the optimizer still regard a solution where 4x + y + z = 10 as a valid solution?

16

Unfortunately this error stops execution during the optimize!() call, so I am unable to check the status as the function does not terminate naturally

17

Unfortunately this error stops execution during the optimize!() call, so I am unable to check the status as the function does not terminate naturally

optimize!(model) should be fine. The error will occur when you call value(x). It probably means that your system has no feasible solution.

18

If you would just like to give a range per variable (rather than a description of all feasible points), you can take any variable, first minimize this variable as the objective, then maximize it. This will give you an interval for each variable. Note that some problems may turn out unbounded, in which case there is no lower/upper bound.

This procedure is also known as OBBT (Optimization-based bound tightening) in the literature.

1 Like
19

In HPolyhedron, strict inequalities are interpreted as non-strict inequalities <=.