Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Is the redistribution available to solve a nonlinear (mixed) complementarity problem? #87

Closed
a24hori opened this issue Jul 21, 2023 · 4 comments · Fixed by #82
Closed

Comments

@a24hori
Copy link

a24hori commented Jul 21, 2023

The original PATH solver enables to solve a nonlinear CP by running pathmcp(...).

@odow
Copy link
Collaborator

odow commented Jul 21, 2023

You can solve nonlinear MCP if you provide the function and Jacobian directly:

PATHSolver.jl/src/C_API.jl

Lines 604 to 752 in 87f8884

"""
solve_mcp(
F::Function,
J::Function
lb::Vector{Cdouble},
ub::Vector{Cdouble},
z::Vector{Cdouble};
nnz::Int = length(lb)^2,
variable_name::Vector{String}=String[],
constraint_name::Vector{String}=String[],
silent::Bool = false,
generate_output::Integer = 0,
use_start::Bool = true,
use_basics::Bool = false,
jacobian_structure_constant::Bool = false,
jacobian_data_contiguous::Bool = false,
jacobian_linear_elements::Vector{Int} = Int[],
kwargs...
)
Mathematically, the mixed complementarity problem is to find an x such that
for each i, at least one of the following hold:
1. F_i(x) = 0, lb_i <= (x)_i <= ub_i
2. F_i(x) > 0, (x)_i = lb_i
3. F_i(x) < 0, (x)_i = ub_i
where F is a given function from R^n to R^n, and lb and ub are prescribed
lower and upper bounds.
## The `F` argument
`F` is a function that calculates the value of function ``F(x)`` and stores the
result in `f`. It must have the signature:
```julia
function F(n::Cint, x::Vector{Cdouble}, f::Vector{Cdouble})
for i in 1:n
f[i] = ... do stuff ...
end
return Cint(0)
end
```
## The `J` argument
`J` is a function that calculates the Jacobiann of the function ``F(x)``. The
Jacobian is a square sparse matrix. It must have the signature:
```julia
function J(
n::Cint,
nnz::Cint,
x::Vector{Cdouble},
col::Vector{Cint},
len::Vector{Cint},
row::Vector{Cint},
data::Vector{Cdouble},
)
# ...
return Cint(0)
end
```
where:
* `n` is the number of variables (which is also the number rows and columns in
the Jacobian matrix).
* `nnz` is the maximum number of non-zero terms in the Jacobian. This value is
chosen by the user as the `nnz` argumennt to `solve_mcp`.
* `x` is the value of the decision variables at which to evaluate the Jacobian.
The remaining arguments, `col`, `len`, `row`, and `data`, specify a sparse
column representation of the Jacobian matrix. These must be filled in by your
function.
* `col` is a length `n` vector, where `col[i]` is the 1-indexed position of the
start of the non-zeros in column `i` in the `data` vector.
* `len` is a length `n` vector, where `len[i]` is the number of non-zeros in
column `i`.
Together, `col` and `len` can be used to form a range of indices in `row` and
`data` corresponding to the non-zero elements of column `i` in the Jacobian.
Thus, we can iterate over the non-zeros in the Jacobian using:
```julia
for i in 1:n
for k in (col[i]):(col[i] + len[i] - 1)
row[k] = ... the 1-indexed row of the k'th non-zero in the Jacobian
data[k] = ... the value of the k'th non-zero in the Jacobian
end
end
```
To improve performance, see the `jacobian_structure_constant` and
`jacobian_data_contiguous` keyword arguments.
## Other positional arguments
* `lb`: a vector of the variable lower bounds
* `ub`: a vector of the variable upper bounds
* `z`: an initial starting point for the search. You can disable this by
passing an empty vector and settig `use_start = false`
## Keyword arguments
* `nnz`: the maximum number of non-zeros in the Jacobian matrix. If not
specified if defaults to the dense estimate of `n^2` where `n` is the number
of variables.
* `variable_name`: a vector of variable names. This can improve the legibility
of the output printed by PATH, particularly if there are issues associated
with a particular variable.
* `constraint_name`: a vector of constraint names. This can improve the
legibility of the output printed by PATH, particularly if there are issues
associated with a particular row of the `F` function.
* `silent`: set `silent = true` to disable printing.
* `generate_output`: an integer mask passed to the C API of PATH to dictate
with output can be displayed.
* `use_start`: set `use_start = false` to disable the use of the startint point
`z`.
* `use_basics`: set `use_basics = true` to use the basis provided.
* `jacobian_structure_constant`: if `true`, the sparsity pattern of the
Jacobian matrix must be constant between evaluations. You can improve
performance by setting this to `true` and filling the `col`, `len` and `row`
on the first evaluation only.
* `jacobian_data_contiguous`: if `true`, the Jacobian data is stored
contiguously from `1..nnz` in the `row` and `data` arrays of the Jacobian
callback. In most cases, you can improve performance by settinng this to
`true`. It is `false` by default for the general case.
* `jacobian_linear_elements`: a vector of the 1-indexed indices of the Jacobian
`data` array that appear linearly in the Jacobian, that is, their value is
independent of the point `x` at which the Jacobian is evaluated. If you set
this option, you must also set `jacobian_structure_constant = true`.
* `kwargs`: other options passed to directly to PATH.
"""
function solve_mcp(
F::Function,
J::Function,
lb::Vector{Cdouble},
ub::Vector{Cdouble},
z::Vector{Cdouble};
nnz::Int = length(lb)^2,
variable_names::Vector{String} = String[],
constraint_names::Vector{String} = String[],
silent::Bool = false,
generate_output::Integer = 0,
use_start::Bool = true,
use_basics::Bool = false,
jacobian_structure_constant::Bool = false,
jacobian_data_contiguous::Bool = false,
jacobian_linear_elements::Vector{Int} = Int[],
kwargs...,
)

An upcoming version of JuMP will support nonlinear complementarity problems: #82

@a24hori
Copy link
Author

a24hori commented Jul 21, 2023

I appreciate your reply. I am a beginner of julia and JuMP, so I am confused some expressions in JuMP. I have one important question to use JuMP.

In the example shown in README, you demonstrate solving the LCP Mx+q \perp x>= 0; however, the standard form of the LCP also has the nonnegative constraint of Mx+q, i.e., Mx+q>=0. Is this automatically added in the model?

@odow
Copy link
Collaborator

odow commented Jul 21, 2023

A complementarity problem in JuMP depends only on the bounds of the x variable:

https://jump.dev/JuMP.jl/stable/moi/reference/standard_form/#MathOptInterface.Complements

For Mx+q ⟂ x>= 0, If x = 0 then Mx+q >= 0, and if x > 0, then Mx+q == 0.

JuMP is not like Pyomo or AMPL, where you can put inequalities on the function as well as the variable.

@hurak
Copy link

hurak commented Aug 2, 2023

My guess is that @a24hori is just confused by reading the complementarity constraint in the form Mx+q ⟂ x>= 0, while the form used elsewhere is 0<=Mx+q ⟂ x>= 0. If I understand it correctly, then indeed, the nonnegativity condition on Mx+q is automatically introduced anyway.

@odow odow closed this as completed in #82 Aug 16, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants