[Bug] Unstable solve on GPU

[jbcaillau]

The following solve on GPU is not stable: successive runs with fixed max_iter may produce different iterates. no such issue on CPU.

NB. problem is probably not feasible. dynamics with x(0) == [1, 2, 3, 4]; (∂(x₁))(t) == x₁(t) ^ 2 + x₂(t) ^ 2 most likely leads to explosion.

julia> o = @def begin
                   t ∈ [0, 1], time
                   x ∈ R⁴, state
                   u ∈ R², control

                   x(0) == [1, 2, 3, 4]

                   ∂(x₁)(t) == x₁(t)^2 + x₂(t)^2
                   ∂(x₂)(t) == u₁(t) * u₂(t)
                   ∂(x₃)(t) == x₂(t) + x₃(t) + x₄(t)
                   ∂(x₄)(t) == u₁(t)

                   (x₁(1) + x₂(1) + x₃(1))^2 + 0.5∫( u₁(t)^2 + u₂(t)^2 ) → min
               end
Abstract definition:

    t ∈ [0, 1], time
    x ∈ R⁴, state
    u ∈ R², control
    x(0) == [1, 2, 3, 4]
    (∂(x₁))(t) == x₁(t) ^ 2 + x₂(t) ^ 2
    (∂(x₂))(t) == u₁(t) * u₂(t)
    (∂(x₃))(t) == x₂(t) + x₃(t) + x₄(t)
    (∂(x₄))(t) == u₁(t)
    (x₁(1) + x₂(1) + x₃(1)) ^ 2 + 0.5 * ∫(u₁(t) ^ 2 + u₂(t) ^ 2) → min

The (autonomous) optimal control problem is of the form:

    minimize  J(x, u) = g(x(0), x(1)) + ∫ f⁰(x(t), u(t)) dt, over [0, 1]

    subject to

        ẋ(t) = f(x(t), u(t)), t in [0, 1] a.e.,

        ϕ₋ ≤ ϕ(x(0), x(1)) ≤ ϕ₊, 

    where x(t) ∈ R⁴ and u(t) ∈ R².


julia> sol = solve(o, :exa, :madnlp; max_iter=10, exa_backend=backend)
▫ This is OptimalControl version v1.1.6 running with: direct, exa, madnlp.

▫ The optimal control problem is solved with CTDirect version v0.17.4.

   ┌─ The NLP is modelled with ExaModels and solved with MadNLP suite.
   │
   ├─ Number of time steps⋅: 250
   └─ Discretisation scheme: midpoint

▫ This is MadNLP version v0.8.12, running with cuDSS v0.7.1

Number of nonzeros in constraint Jacobian............:     4254
Number of nonzeros in Lagrangian Hessian.............:     5006

Total number of variables............................:     1506
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1004
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du inf_compl lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.0000000e-01 3.90e+00 6.00e-01 1.58e-07  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.0000355e-01 3.90e+00 8.28e-03 3.48e-06  -1.0 1.88e+02    -  1.00e+00 2.00e-07h  2
   2  9.7919791e+00 2.92e+00 8.88e-04 3.11e-06  -1.0 2.70e+01    -  1.00e+00 2.50e-01h  3
   3  5.3284602e+01 1.60e-01 5.17e-04 2.00e-06  -1.0 9.50e+00   0.0 1.00e+00 1.00e+00h  1
   4  1.0994023e+02 2.45e-02 6.47e-03 2.03e-06  -2.5 3.30e+00   0.4 1.00e+00 1.00e+00h  1
   5  1.7675944e+02 2.19e-02 6.29e-02 2.23e-06  -3.8 2.65e+00  -0.1 1.00e+00 1.00e+00h  1
   6  2.0476903e+02 1.04e-01 3.21e+01 2.41e-06  -3.8 1.89e+01  -0.5 1.00e+00 2.50e-01h  3
   7  2.5090064e+02 8.95e-02 9.29e+01 2.81e-06  -3.8 6.93e+00  -0.1 1.00e+00 2.50e-01h  3
   8  3.1185887e+02 1.00e-01 1.00e+02 2.93e-06  -3.8 1.16e+01  -0.6 1.00e+00 2.50e-01h  3
   9  4.0594684e+02 8.04e-02 2.50e+02 3.63e-06  -3.8 1.06e+01  -0.2 1.00e+00 2.38e-01h  3
iter    objective    inf_pr   inf_du inf_compl lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  4.7149148e+02 7.24e-02 3.51e+02 4.29e-06  -3.8 1.38e+01  -0.6 1.00e+00 1.25e-01h  4

Number of Iterations....: 10

                                   (scaled)                 (unscaled)
Objective...............:   4.7149148065749313e+02    4.7149148065749313e+02
Dual infeasibility......:   3.5055929454207512e+02    3.5055929454207512e+02
Constraint violation....:   7.2358492745758812e-02    7.2358492745758812e-02
Complementarity.........:   4.2945788703848589e-06    4.2945788703848589e-06
Overall NLP error.......:   3.5055929454207512e+02    3.5055929454207512e+02

Number of objective function evaluations             = 31
Number of objective gradient evaluations             = 11
Number of constraint evaluations                     = 31
Number of constraint Jacobian evaluations            = 11
Number of Lagrangian Hessian evaluations             = 10
Total wall-clock secs in solver (w/o fun. eval./lin. alg.)  =  0.083
Total wall-clock secs in linear solver                      =  0.048
Total wall-clock secs in NLP function evaluations           =  0.012
Total wall-clock secs                                       =  0.143

EXIT: Maximum Number of Iterations Exceeded.
• Solver:
  ✓ Successful  : false
  │  Status     : MAXIMUM_ITERATIONS_EXCEEDED
  │  Message    : MadNLP
  │  Iterations : 10
  │  Objective  : 471.4914806574931
  └─ Constraints violation : 0.07235849274575881

• Boundary duals: [0.0, 0.0, 0.0, 0.0]


julia> sol = solve(o, :exa, :madnlp; max_iter=10, exa_backend=backend)
▫ This is OptimalControl version v1.1.6 running with: direct, exa, madnlp.

▫ The optimal control problem is solved with CTDirect version v0.17.4.

   ┌─ The NLP is modelled with ExaModels and solved with MadNLP suite.
   │
   ├─ Number of time steps⋅: 250
   └─ Discretisation scheme: midpoint

▫ This is MadNLP version v0.8.12, running with cuDSS v0.7.1

Number of nonzeros in constraint Jacobian............:     4254
Number of nonzeros in Lagrangian Hessian.............:     5006

Total number of variables............................:     1506
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1004
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du inf_compl lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.0000000e-01 3.90e+00 6.00e-01 1.58e-07  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.0000355e-01 3.90e+00 8.28e-03 3.48e-06  -1.0 1.88e+02    -  1.00e+00 2.00e-07h  2
   2  9.7919769e+00 2.92e+00 8.88e-04 3.11e-06  -1.0 2.70e+01    -  1.00e+00 2.50e-01h  3
   3  5.3284532e+01 1.60e-01 5.17e-04 2.00e-06  -1.0 9.50e+00   0.0 1.00e+00 1.00e+00h  1
   4  1.0994010e+02 2.45e-02 6.47e-03 2.03e-06  -2.5 3.30e+00   0.4 1.00e+00 1.00e+00h  1
   5  1.7675926e+02 2.19e-02 6.29e-02 2.23e-06  -3.8 2.65e+00  -0.1 1.00e+00 1.00e+00h  1
   6  2.0476884e+02 1.04e-01 3.21e+01 2.41e-06  -3.8 1.89e+01  -0.5 1.00e+00 2.50e-01h  3
   7  2.5090041e+02 8.95e-02 9.29e+01 2.81e-06  -3.8 6.93e+00  -0.1 1.00e+00 2.50e-01h  3
   8  3.1185857e+02 1.00e-01 1.00e+02 2.93e-06  -3.8 1.16e+01  -0.6 1.00e+00 2.50e-01h  3
   9  4.0594652e+02 8.04e-02 2.50e+02 3.63e-06  -3.8 1.06e+01  -0.2 1.00e+00 2.38e-01h  3
iter    objective    inf_pr   inf_du inf_compl lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  4.7149111e+02 7.24e-02 3.51e+02 4.29e-06  -3.8 1.38e+01  -0.6 1.00e+00 1.25e-01h  4

Number of Iterations....: 10

                                   (scaled)                 (unscaled)
Objective...............:   4.7149110750267818e+02    4.7149110750267818e+02
Dual infeasibility......:   3.5055899993134869e+02    3.5055899993134869e+02
Constraint violation....:   7.2358454841026612e-02    7.2358454841026612e-02
Complementarity.........:   4.2945774528465284e-06    4.2945774528465284e-06
Overall NLP error.......:   3.5055899993134869e+02    3.5055899993134869e+02

Number of objective function evaluations             = 31
Number of objective gradient evaluations             = 11
Number of constraint evaluations                     = 31
Number of constraint Jacobian evaluations            = 11
Number of Lagrangian Hessian evaluations             = 10
Total wall-clock secs in solver (w/o fun. eval./lin. alg.)  =  0.209
Total wall-clock secs in linear solver                      =  0.044
Total wall-clock secs in NLP function evaluations           =  0.012
Total wall-clock secs                                       =  0.264

EXIT: Maximum Number of Iterations Exceeded.
• Solver:
  ✓ Successful  : false
  │  Status     : MAXIMUM_ITERATIONS_EXCEEDED
  │  Message    : MadNLP
  │  Iterations : 10
  │  Objective  : 471.4911075026782
  └─ Constraints violation : 0.07235845484102661

• Boundary duals: [0.0, 0.0, 0.0, 0.0]


julia> sol = solve(o, :exa, :madnlp; max_iter=10, exa_backend=backend)
▫ This is OptimalControl version v1.1.6 running with: direct, exa, madnlp.

▫ The optimal control problem is solved with CTDirect version v0.17.4.

   ┌─ The NLP is modelled with ExaModels and solved with MadNLP suite.
   │
   ├─ Number of time steps⋅: 250
   └─ Discretisation scheme: midpoint

▫ This is MadNLP version v0.8.12, running with cuDSS v0.7.1

Number of nonzeros in constraint Jacobian............:     4254
Number of nonzeros in Lagrangian Hessian.............:     5006

Total number of variables............................:     1506
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1004
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du inf_compl lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.0000000e-01 3.90e+00 6.00e-01 1.58e-07  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.0000355e-01 3.90e+00 8.28e-03 3.48e-06  -1.0 1.88e+02    -  1.00e+00 2.00e-07h  2
   2  9.8103812e+00 2.92e+00 9.08e-04 3.11e-06  -1.0 2.68e+01  -4.0 1.00e+00 2.50e-01h  3
   3  2.2668124e+01 1.46e+00 3.75e-03 2.56e-06  -1.0 1.02e+01  -0.9 1.00e+00 5.00e-01h  2
   4  1.2942906e+02 1.09e-01 7.33e-05 2.00e-06  -1.0 5.85e+00  -0.4 1.00e+00 1.00e+00h  1
   5  1.5004709e+02 1.84e-01 3.39e-04 2.01e-06  -2.5 6.99e+00  -0.9 1.00e+00 1.00e+00h  1
   6  1.9524911e+02 9.50e-02 6.84e-01 2.01e-06  -2.5 2.65e+00  -0.5 1.00e+00 5.00e-01h  2
   7  2.0983066e+02 8.62e-02 2.59e+00 2.02e-06  -2.5 7.51e+00  -1.0 1.00e+00 1.25e-01h  4
   8  2.5814608e+02 6.90e-02 4.01e+00 2.03e-06  -2.5 6.95e+00  -0.5 1.00e+00 2.50e-01h  3
   9  2.9647696e+02 6.08e-02 1.28e+01 2.09e-06  -2.5 9.31e+00  -0.1 1.00e+00 1.25e-01h  4
iter    objective    inf_pr   inf_du inf_compl lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  3.3940577e+02 5.46e-02 1.19e+01 2.09e-06  -2.5 1.06e+01  -0.6 1.00e+00 1.25e-01h  4

Number of Iterations....: 10

                                   (scaled)                 (unscaled)
Objective...............:   3.3940577196717010e+02    3.3940577196717010e+02
Dual infeasibility......:   1.1944768937354576e+01    1.1944768937354576e+01
Constraint violation....:   5.4637889678052810e-02    5.4637889678052810e-02
Complementarity.........:   2.0883827423881080e-06    2.0883827423881080e-06
Overall NLP error.......:   1.1944768937354576e+01    1.1944768937354576e+01

Number of objective function evaluations             = 34
Number of objective gradient evaluations             = 11
Number of constraint evaluations                     = 34
Number of constraint Jacobian evaluations            = 11
Number of Lagrangian Hessian evaluations             = 10
Total wall-clock secs in solver (w/o fun. eval./lin. alg.)  =  0.094
Total wall-clock secs in linear solver                      =  0.064
Total wall-clock secs in NLP function evaluations           =  0.014
Total wall-clock secs                                       =  0.172

EXIT: Maximum Number of Iterations Exceeded.
• Solver:
  ✓ Successful  : false
  │  Status     : MAXIMUM_ITERATIONS_EXCEEDED
  │  Message    : MadNLP
  │  Iterations : 10
  │  Objective  : 339.4057719671701
  └─ Constraints violation : 0.05463788967805281

• Boundary duals: [0.0, 0.0, 0.0, 0.0]


julia> sol = solve(o, :exa, :madnlp; max_iter=10, exa_backend=backend)
▫ This is OptimalControl version v1.1.6 running with: direct, exa, madnlp.

▫ The optimal control problem is solved with CTDirect version v0.17.4.

   ┌─ The NLP is modelled with ExaModels and solved with MadNLP suite.
   │
   ├─ Number of time steps⋅: 250
   └─ Discretisation scheme: midpoint

▫ This is MadNLP version v0.8.12, running with cuDSS v0.7.1

Number of nonzeros in constraint Jacobian............:     4254
Number of nonzeros in Lagrangian Hessian.............:     5006

Total number of variables............................:     1506
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1004
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du inf_compl lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.0000000e-01 3.90e+00 6.00e-01 1.58e-07  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.0000355e-01 3.90e+00 8.28e-03 3.48e-06  -1.0 1.88e+02    -  1.00e+00 2.00e-07h  2
   2  5.1449485e+01 1.95e+00 1.49e-03 2.74e-06  -1.0 1.82e+01  -2.0 1.00e+00 5.00e-01h  2
   3  7.9293991e+01 5.24e-01 3.21e-05 2.00e-06  -1.0 1.19e+01  -0.7 1.00e+00 1.00e+00h  1
   4  1.2006318e+02 1.19e-02 2.89e-04 2.00e-06  -1.7 3.52e+00  -0.2 1.00e+00 1.00e+00h  1
   5  1.5471721e+02 1.63e-02 9.18e+00 2.17e-06  -3.8 5.77e+00  -0.7 1.00e+00 5.00e-01h  2
   6  2.2239687e+02 3.32e-02 3.09e+01 2.23e-06  -3.8 6.92e+00  -0.3 1.00e+00 5.00e-01h  2
   7  2.7195615e+02 5.32e-02 7.82e+01 2.60e-06  -3.8 1.13e+01  -0.8 1.00e+00 2.50e-01h  3
   8  3.1523638e+02 4.85e-02 2.32e+02 3.77e-06  -3.8 1.06e+01  -0.3 1.00e+00 1.25e-01h  4
   9  3.6219801e+02 4.81e-02 2.29e+02 3.92e-06  -3.8 2.06e+01  -0.8 1.00e+00 1.25e-01h  4
iter    objective    inf_pr   inf_du inf_compl lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  4.2547402e+02 4.75e-02 3.73e+02 4.80e-06  -3.8 1.35e+01  -0.4 1.00e+00 1.25e-01h  4

Number of Iterations....: 10

                                   (scaled)                 (unscaled)
Objective...............:   4.2547401691846102e+02    4.2547401691846102e+02
Dual infeasibility......:   3.7297322993663613e+02    3.7297322993663613e+02
Constraint violation....:   4.7491506630344779e-02    4.7491506630344779e-02
Complementarity.........:   4.7952429881045119e-06    4.7952429881045119e-06
Overall NLP error.......:   3.7297322993663613e+02    3.7297322993663613e+02

Number of objective function evaluations             = 33
Number of objective gradient evaluations             = 11
Number of constraint evaluations                     = 33
Number of constraint Jacobian evaluations            = 11
Number of Lagrangian Hessian evaluations             = 10
Total wall-clock secs in solver (w/o fun. eval./lin. alg.)  =  0.094
Total wall-clock secs in linear solver                      =  0.072
Total wall-clock secs in NLP function evaluations           =  0.014
Total wall-clock secs                                       =  0.179

EXIT: Maximum Number of Iterations Exceeded.
• Solver:
  ✓ Successful  : false
  │  Status     : MAXIMUM_ITERATIONS_EXCEEDED
  │  Message    : MadNLP
  │  Iterations : 10
  │  Objective  : 425.474016918461
  └─ Constraints violation : 0.04749150663034478

• Boundary duals: [0.0, 0.0, 0.0, 0.0]

[jbcaillau]

@amontoison any clue?

[jbcaillau]

• @jbcaillau re-check with CTParser fix (scalar vs. vector of length one bounds on exa)

EDIT: the fix solves indeed random NLP constraint numbererror (both on CPU and GPU), but sill unstable runs on GPU (minor changes, but changes on iterations number and then objective, etc.)

[PierreMartinon]

@jbcaillau maybe this is due to the problem and settings modifications on GPU, cf https://madsuite.org/MadNLP.jl/dev/tutorials/gpu/
Can you try to set explicitely tol and bound_relax_factor to a lower value eg 1e-7 as suggested in the example ?

[jbcaillau]

cc: @Dtinas10

[jbcaillau]

@PierreMartinon indeed, stable solve (on an example stemming from quantum control) with the following call:

sol = solve(o, :exa, :madnlp; init=(control=φ_init,), tol=tol, grid_size=N, exa_backend=CUDABackend(), bound_relax_factor=tol)  # GPU

still don't see what can change from one call to the other without this option as everything should be deterministic (as on CPU where runs are absolutely stable)

@amontoison any clue?

@jbcaillau maybe this is due to the problem and settings modifications on GPU, cf https://madsuite.org/MadNLP.jl/dev/tutorials/gpu/ Can you try to set explicitely tol and bound_relax_factor to a lower value eg 1e-7 as suggested in the example ?

[jbcaillau]

@PierreMartinon NB: N = 10000 on A100 🔥

julia> sol = solve(o, :exa, :madnlp; init=(control=φ_init,), tol=tol, grid_size=N, exa_backend=CUDABackend(), bound_relax_factor=tol)  # GPU
▫ This is OptimalControl version v1.1.8-beta.3 running with: direct, exa, madnlp.

▫ The optimal control problem is solved with CTDirect version v0.17.5-beta.

   ┌─ The NLP is modelled with ExaModels and solved with MadNLP suite.
   │
   ├─ Number of time steps⋅: 10000
   └─ Discretisation scheme: midpoint

▫ This is MadNLP version v0.8.12, running with cuDSS v0.7.1

Number of nonzeros in constraint Jacobian............:  2880032
Number of nonzeros in Lagrangian Hessian.............:  4480064

Total number of variables............................:   330033
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:   320032
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du inf_compl lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  9.9000000e-01 9.00e-01 5.00e-02 3.96e-07  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  9.9010072e-01 9.00e-01 2.71e-02 3.96e-05  -1.7 3.45e+01  -4.0 1.00e+00 2.00e-05h  1
   2 -6.5303999e-02 1.23e-04 7.36e-08 2.00e-05  -1.7 1.12e+00  -4.5 1.00e+00 1.00e+00h  1
   3 -7.2520473e-03 3.61e-05 2.39e-05 3.11e-06  -5.7 6.89e-01  -5.0 1.00e+00 1.00e+00h  1
   4  9.6570644e-04 5.24e-05 2.49e-05 1.86e-06  -5.7 9.44e-01    -  1.00e+00 1.00e+00h  1
   5  8.4842599e-04 3.04e-06 1.53e-06 1.84e-06  -5.7 1.97e-01  -5.4 1.00e+00 1.00e+00h  1
   6  4.1752310e-04 3.07e-06 3.12e-04 9.46e-08  -8.0 1.98e-01    -  9.79e-01 1.00e+00h  1
   7 -3.4017124e-03 3.84e-07 5.29e-04 2.35e-08  -8.0 7.45e-02  -5.9 8.96e-01 1.00e+00h  1
   8 -4.2961675e-03 6.71e-07 3.39e-07 1.14e-08  -8.0 9.59e-02    -  1.00e+00 1.00e+00h  1
   9 -4.2249737e-03 1.88e-08 9.60e-09 9.15e-09  -8.0 1.71e-02  -6.4 1.00e+00 1.00e+00h  1

Number of Iterations....: 9

                                   (scaled)                 (unscaled)
Objective...............:  -4.2249736522177450e-03   -4.2249736522177450e-03
Dual infeasibility......:   9.6029412527972555e-09    9.6029412527972555e-09
Constraint violation....:   1.8816179336913941e-08    1.8816179336913941e-08
Complementarity.........:   9.1499574525742553e-09    9.1499574525742553e-09
Overall NLP error.......:   1.8816179336913941e-08    1.8816179336913941e-08

Number of objective function evaluations             = 10
Number of objective gradient evaluations             = 10
Number of constraint evaluations                     = 10
Number of constraint Jacobian evaluations            = 10
Number of Lagrangian Hessian evaluations             = 9
Total wall-clock secs in solver (w/o fun. eval./lin. alg.)  =  0.798
Total wall-clock secs in linear solver                      =  0.092
Total wall-clock secs in NLP function evaluations           =  0.035
Total wall-clock secs                                       =  0.925

EXIT: Optimal Solution Found (tol = 1.0e-07).
• Solver:
  ✓ Successful  : true
  │  Status     : SOLVE_SUCCEEDED
  │  Message    : MadNLP
  │  Iterations : 9
  │  Objective  : -0.004224973652217745
  └─ Constraints violation : 1.881617933691394e-8

• Boundary duals: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

[amontoison]

@PierreMartinon indeed, stable solve (on an example stemming from quantum control) with the following call:

sol = solve(o, :exa, :madnlp; init=(control=φ_init,), tol=tol, grid_size=N, exa_backend=CUDABackend(), bound_relax_factor=tol) # GPU
still don't see what can change from one call to the other without this option as everything should be deterministic (as on CPU where runs are absolutely stable)

@amontoison any clue?

@jbcaillau maybe this is due to the problem and settings modifications on GPU, cf https://madsuite.org/MadNLP.jl/dev/tutorials/gpu/ Can you try to set explicitely tol and bound_relax_factor to a lower value eg 1e-7 as suggested in the example ?

@jbcaillau with bound_relax_factor=tol we relax all equality g(x) = 0 into inequality constraints -tol <= g(x) <= tol.
By default we solve with an accuracy of 1e-4 on GPU and the default value for the option bound_relax_factor is a such that the relaxation of the constraints will not really impact the problem solution for the accuracy requested.
BUT if you want a more accurate solution, than you need to adapt the value of bound_relax_factor.

[jbcaillau]

@amontoison the relaxation process is clear, thanks. the question on top of that is why would there be some unstable = stochastic behaviour, that is successive runs on GPU giving slightly different iterations and results (this is not happening on CPU). setting bound_relax_factor seems to fix that on the particular example i mentioned, but that’s does not explain anything.

@jbcaillau with bound_relax_factor=tol we relax all equality g(x) = 0 into inequality constraints -tol <= g(x) <= tol. By default we solve with an accuracy of 1e-4 on GPU and the default value for the option bound_relax_factor is a such that the relaxation of the constraints will not really impact the problem solution for the accuracy requested. BUT if you want a more accurate solution, than you need to adapt the value of bound_relax_factor.

[jbcaillau]

@PierreMartinon we also need to check https://madsuite.org/MadNLP.jl/dev/tutorials/gpu/#Solving-the-problem-on-the-GPU-with-HyKKT

[amontoison]

@jbcaillau CUDSS.jl is not deterministic (by default) and more the linear system is ill-conditioned, more risks you have to get different results between two runs.
The option to relax the equality constraints could help the stability of KKT system and make it less sensible to non-determinism.

[amontoison]

If you want a stable solver on GPU, MadNCL.jl is what you need 🙂
The plan to use it in the soon future at Los Alamos instead of MadNLP.jl for gaz problems (that have nasty complementarity constraints).

[amontoison]

If you want a stable solver on GPU, MadNCL.jl is what you need 🙂
The plan to use it in the soon future at Los Alamos instead of MadNLP.jl for gaz problems (that have nasty complementarity constraints).