jump-dev / Clp.jl

A Julia interface to the Coin-OR Linear Programming solver (CLP)
https://projects.coin-or.org/Clp
Other
54 stars 25 forks source link

cannot get the primal solution when using the Clp ipm method #125

Closed nht2018 closed 2 years ago

nht2018 commented 2 years ago

Hello, I came across two problems.

1

After I use iterior point method(ipm) to solve the LP programming problem, I cannot get the primal solution through value.(x). (But there is no problem with simplex method.) Here are my codes:

using JuMP
using Clp

#read data
...

model = Model(Clp.Optimizer);
@variable(model, x[1:n]);
for i in 1: n
    set_lower_bound(x[i], lb[i]);
    set_upper_bound(x[i], ub[i]);
end
@constraint(model, constraint, A * x .== b);
@objective(model, Min, dot(x, c));
set_optimizer_attribute(model, "SolveType", 4);
set_optimizer_attribute(model, "PresolveType", 1);   
set_optimizer_attribute(model, "MaximumSeconds", 1e3);
set_optimizer_attribute(model, "PrimalTolerance", 1e-6);
set_optimizer_attribute(model, "DualTolerance", 1e-6);
optimize!(model);

where A is a m*n matrix, c, lb, ub are n*1 matrices, b is m*1 matrix. When I try to get the primal solution by value.(x), it traces back:

ERROR: Primal solution not available

2

I have tried different approaches to get the iteration counts, including

MOI.get(model, MOI.SimplexItertions()) 
solution_summary(model).simplex_iterations
solution_summary(model).barrier_iterations
Clp_numberIterations(model)

but none works.

Could you please provide some approaches to acquire the iteration counts and primal solution? Thanks very much.

odow commented 2 years ago

Did it find a feasible solution?

After optimize!(model) you must check termination_status(model) and primal_status(model). What are they?

Please provide a reproducible example.

nht2018 commented 2 years ago

Thanks for your quick reply. The termination_status is OPTIMAL and I could get the dual solution by using dual.(constraint). And primal_status is UNKNOWN_RESULT_STATUS.

The attachmments are my codes, one instance(cont1 in the dataset mittelmann, which has been prosolved), and the log printed out.(It seems that the algorithm works well.)

example.zip

odow commented 2 years ago

If the primal status is UNKNOWN_RESULT_STATUS then something went wrong in Clp. For some reason, Clp didn't declare Cbc_primalFeasible: https://github.com/jump-dev/Clp.jl/blob/eb6ec27e1207593d98671dbdc2b429d3f2961d47/src/MOI_wrapper/MOI_wrapper.jl#L437-L440

This is not a problem with Clp.jl. But you'd have to dig around in Clp to understand what happened in more details. Can Clp solve this via an MPS file? Is it a tolerance issue? It it your presolve routine missing with scaling?

nht2018 commented 2 years ago

I find the problem. It is necessary to take crossover steps to get the primal solution. After I change

set_optimizer_attribute(model, "SolveType", 4);  # ipm without crossover

to

set_optimizer_attribute(model, "SolveType", 3);  # ipm with crossover

then value.(x) returns the primal solution.

By the way, I have another question: how to get the iteration counts? I have tried

MOI.get(model, MOI.SimplexItertions()) 
solution_summary(model).simplex_iterations
solution_summary(model).barrier_iterations
Clp_numberIterations(model)

but none of them works.

odow commented 2 years ago

Potentially:

model = Model(Clp.Optimizer)
...
clp = unsafe_backend(model)
Clp_numberIterations(clp)

but that might just be for simplex.

If it doesn't work, there might be something else in the C API: https://github.com/jump-dev/Clp.jl/blob/master/src/libClp.jl

nht2018 commented 2 years ago

It works for both simplex and ipm method. Thanks very much.

odow commented 2 years ago

Great. I'll close this as fixed. If you have further questions with your benchmarking, please post on the community forum: https://discourse.julialang.org/c/domain/opt/13