diff --git a/.gitignore b/.gitignore index fc3a525f..2cc34922 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,5 @@ Manifest.toml examples/.ipynb_checkpoints .DS_Store /Manifest.toml +*.sbatch +*.out diff --git a/Project.toml b/Project.toml index d73a81b4..21862f2b 100644 --- a/Project.toml +++ b/Project.toml @@ -5,13 +5,26 @@ version = "0.4.2" [deps] BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0" +CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" +ClusterManagers = "34f1f09b-3a8b-5176-ab39-66d58a4d544e" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" MathOptSetDistances = "3b969827-a86c-476c-9527-bb6f1a8fbad5" +NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd" +PGLib = "07a8691f-3d11-4330-951b-3c50f98338be" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +PowerModels = "c36e90e8-916a-50a6-bd94-075b64ef4655" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] diff --git a/bilevel.jl b/bilevel.jl new file mode 100644 index 00000000..75dd30fe --- /dev/null +++ b/bilevel.jl @@ -0,0 +1,131 @@ +using Distributed + +# Add worker processes if needed +# addprocs() + +@everywhere pkg_path = @__DIR__ + +@everywhere import Pkg + +@everywhere Pkg.activate(pkg_path) + +@everywhere Pkg.instantiate() + +@everywhere using BilevelJuMP +@everywhere using Ipopt +@everywhere using QuadraticToBinary +@everywhere using Gurobi +@everywhere using LinearAlgebra +@everywhere using Combinatorics +@everywhere using StatsBase +@everywhere using Random + +@everywhere function build_lower_level(num_nodes=10; jump_model=Model(), seed=1, qS=100) + Random.seed!(seed) + costs = rand(1:1:200, num_nodes) + demands = rand(1:1:100, num_nodes) + # circle network with num_nodes extra random lines + comb = collect(combinations(1:num_nodes,2)) + num_comb = length(comb) + rand_lines = comb[sample(1:num_comb, ceil(Int, num_nodes / 2), replace = false)] + lines = vcat([(i, i+1) for i=1:num_nodes-1], [(num_nodes, 1)], rand_lines) + num_lines = length(lines) + line_limits = rand(1:1:100, num_lines) + @variable(jump_model, 0 <= gS[i=1:num_nodes] <= 100) + @variable(jump_model, 0 <= gR[i=1:num_nodes] <= 40) + @variable(jump_model, 0 <= gD[i=1:num_nodes] <= 100) + # line flows + @variable(jump_model, 0 <= f[i=1:num_lines] <= line_limits[i]) + @objective(jump_model, Min, dot(costs, gR) + 1000 * sum(gD)) + @constraint(jump_model, gS .<= qS) + # demand_equilibrium + # @constraint(jump_model, demand_equilibrium, gS + sum(gR) + gD == 100) + @constraint(jump_model, demand_equilibrium[i=1:num_nodes], + gS[i] + gR[i] + gD[i] + sum(f[j] for j in 1:num_lines if lines[j][1] == i) - sum(f[j] for j in 1:num_lines if lines[j][2] == i) == demands[i] + ) + return jump_model, demand_equilibrium, gS +end + +@everywhere function build_bilevel(num_nodes=10; seed=1) + model = BilevelModel() + @variable(Upper(model), 0 <= qS[i=1:num_nodes] <= 100) + _, demand_equilibrium, gS, = build_lower_level(num_nodes; jump_model=Lower(model), seed=seed, qS=qS) + @variable(Upper(model), lambda[i=1:num_nodes], DualOf(demand_equilibrium[i])) + @objective(Upper(model), Max, dot(lambda, gS)) + return model, qS, lambda +end + +@everywhere function solve_nlp(model) + BilevelJuMP.set_mode(model, BilevelJuMP.StrongDualityMode()) + set_optimizer(model, Ipopt.Optimizer) + start_time = time() + optimize!(model) + end_time = time() + return end_time - start_time +end + +@everywhere function solve_bilevel_exact(model, lambda) + set_optimizer(model, + ()->QuadraticToBinary.Optimizer{Float64}(MOI.instantiate(optimizer_with_attributes(Gurobi.Optimizer, "MIPGap" => 0.1))) + ) + BilevelJuMP.set_mode(model, + BilevelJuMP.FortunyAmatMcCarlMode(dual_big_M = 1000) + ) + set_lower_bound.(lambda, 0.0) + set_upper_bound.(lambda, 1000.0) + start_time = time() + optimize!(model) + end_time = time() + return end_time - start_time +end + +@everywhere function calculate_profit(evaluator_lower_level, demand_equilibrium, gS, _qS, qS) + fix.(_qS, value.(qS)) + # @constraint(evaluator_lower_level, _qS .== value.(qS)) + optimize!(evaluator_lower_level) + return dot(value.(gS), dual.(demand_equilibrium)) +end + +@everywhere function compare_methods(num_nodes=10; seed=1) + # try + bilevel_model, qS, lambda = build_bilevel(num_nodes; seed=seed) + evaluator_lower_level = JuMP.Model(Gurobi.Optimizer) + @variable(evaluator_lower_level, _qS[i=1:num_nodes]) + evaluator_lower_level, demand_equilibrium, gS = build_lower_level(num_nodes; seed=seed, jump_model=evaluator_lower_level, qS=_qS) + nlp_time = solve_nlp(bilevel_model) + opf_time = time() + nlp_profit = calculate_profit(evaluator_lower_level, demand_equilibrium, gS, _qS, qS) + opf_time = time() - opf_time + exact_time = solve_bilevel_exact(bilevel_model, lambda) + # exact_profit = calculate_profit(evaluator_lower_level, demand_equilibrium, gS, _qS, qS) + return nlp_time / opf_time, exact_time / opf_time + # catch e + # return NaN, NaN + # end +end + +nodes = 7:5:57 +results = pmap(x->compare_methods(x), nodes) + +# success nodes +_results = [(i, (r[1], r[2])) for (i,r) in enumerate(results) if !isnan(r[1]) && !isnan(r[2])] +nodes = [nodes[r[1]] for r in _results] +results = [r[2] for r in _results] + +# Save results +using CSV, DataFrames +df = DataFrame(nodes=nodes, nlp=[r[1] for r in results], exact=[r[2] for r in results]) +if isfile(joinpath(@__DIR__, "results", "complexity_mpec.csv")) + df_old = CSV.read(joinpath(@__DIR__, "results", "complexity_mpec.csv"), DataFrame) + df = vcat(df_old, df) +end +CSV.write(joinpath(@__DIR__, "results", "complexity_mpec.csv"), df) + +using Plots +# log y axis scale +plt = Plots.plot(df.nodes, df.nlp, label="NLP", xlabel="Number of nodes", ylabel="Time (x OPF)", yscale=:log10) +Plots.plot!(plt, df.nodes, df.exact, label="Exact", xlabel="Number of nodes", ylabel="Time (x OPF)", yscale=:log10) +Plots.savefig(plt, joinpath(@__DIR__, "results", "complexity_mpec.pdf")) + + + diff --git a/bilevel_diffopt.jl b/bilevel_diffopt.jl new file mode 100644 index 00000000..2747e0fd --- /dev/null +++ b/bilevel_diffopt.jl @@ -0,0 +1,298 @@ +using JuMP +using SparseArrays +using LinearAlgebra +using Ipopt + +################################################ +#= +# Test Linear Bilevel +=# +################################################ + +function constraint_upper!(model, x, y) + @constraints(model, begin + x <= 5 + y <= 8 + y >= 0 + end) +end + +function obj_upper!(model, x, y) + @objective(model, Min, 3x + y) +end + +function build_upper!(model, x, y) + obj_upper!(model, x, y) + constraint_upper!(model, x, y) +end + +function constraint_lower!(model, x, y; slack = zeros(4)) + @constraints(model, begin + x + y + slack[1] <= 8 + 4x + y + slack[2] >= 8 + 2x + y + slack[3] <= 13 + 2x - 7y + slack[4] <= 0 + end) +end + +function obj_lower!(model, x, y, slack) + @objective(model, Min, -x + 10 * sum(s^2 for s in slack)) +end + +function build_lower!(model, x, y, slack) + obj_lower!(model, x, y, slack) + constraint_lower!(model, x, y, slack = slack) +end + +function memoize(foo::Function) + last_x, last_f = nothing, nothing + last_dx, last_dfdx = nothing, nothing + function foo_i(x::T...) where {T<:Real} + if T == Float64 + if x !== last_x + last_x, last_f = x, foo(x...) + end + return last_f::T + else + if x !== last_dx + last_dx, last_dfdx = x, foo(x...) + end + return last_dfdx::T + end + end + return (x...) -> foo_i(x...) +end + +# using BilevelJuMP + +# model = BilevelModel(Ipopt.Optimizer, mode = BilevelJuMP.ProductMode(1e-5)) + +# @variable(Lower(model), x_b) +# @variable(Upper(model), y_b) + +# build_upper!(Upper(model), x_b, y_b) +# build_lower!(Lower(model), x_b, y_b, zeros(4)) +# optimize!(model) + +# objective_value(model) # = 3 * (3.5 * 8/15) + 8/15 # = 6.13... +# value(x_b) # = 3.5 * 8/15 # = 1.86... +# value(y_b) # = 8/15 # = 0.53... + +function test_bilevel_linear() + model_upper = Model(Ipopt.Optimizer) + model_lower = Model(Ipopt.Optimizer) + + @variable(model_upper, y) + @variable(model_upper, x_star) + @variable(model_upper, x_aux) + @variable(model_lower, x) + @variable(model_lower, y_p ∈ MOI.Parameter(1.0)) + @variable(model_lower, slack[1:4]) + + build_upper!(model_upper, x_star, y) + constraint_lower!(model_upper, x_aux, y) + build_lower!(model_lower, x, y_p, slack) + primal_vars = [x; slack] + params = [y_p] + evaluator, cons = create_evaluator(model_lower; x=[primal_vars; params]) + + # Define `f(y) = x_star` & `f'(y)` + function f(y_val) + set_parameter_value(y_p, y_val) + optimize!(model_lower) + @assert is_solved_and_feasible(model_lower) + return value(x) + end + + function ∇f0(y_val) + @assert value(y_p) == y_val + Δs, sp = compute_sensitivity(evaluator, cons, [0.00]; primal_vars=primal_vars, params=params) + return Δs[1] + end + + function ∇f(y_val) + @assert value(y_p) == y_val + Δs, sp = compute_sensitivity(evaluator, cons, [0.001]; primal_vars=primal_vars, params=params) + return Δs[1] + end + + memoized_f = memoize(f) + + @operator(model_upper, op_f, 1, memoized_f, ∇f) + @constraint(model_upper, x_star == op_f(y)) + + optimize!(model_upper) + + @test objective_value(model_upper) ≈ 6.13 atol=0.05 + @test value(x_star) ≈ 1.86 atol=0.05 + @test value(y) ≈ 0.53 atol=0.05 +end + +################################################ +#= +# Test Non-Linear Bilevel +=# +################################################ +function test_bilevel_nonlinear() + upper_model = Model(Ipopt.Optimizer) + lower_model = Model(Ipopt.Optimizer) + + @variable(upper_model, x >= 2) + @variable(lower_model, x_p ∈ MOI.Parameter(1.0)) + + @variable(lower_model, 3 <= y <= 5) + @variable(upper_model, y_star) + + @objective(upper_model, Min, x^4 - sin(y_star)) + + @constraint(upper_model, x^3 + y_star^3 <= 1000) + + @objective(lower_model, Min, y^2 + x_p) + + # Define `f(x_p) = y_star` & `f'(x_p)` + primal_vars = [y] + params = [x_p] + evaluator, cons = create_evaluator(lower_model; x=[primal_vars; params]) + + function f(x_p_val) + set_parameter_value(x_p, x_p_val) + optimize!(lower_model) + @assert is_solved_and_feasible(lower_model) + return value(y) + end + + function ∇f(x_p_val) + @assert value(x_p) == x_p_val + Δs, sp = compute_sensitivity(evaluator, cons, [0.001]; primal_vars=primal_vars, params=params) + return Δs[1] + end + + function ∇f0(x_p_val) + @assert value(x_p) == x_p_val + Δs, sp = compute_sensitivity(evaluator, cons, [0.00]; primal_vars=primal_vars, params=params) + return Δs[1] + end + + memoized_f = memoize(f) + + @operator(upper_model, op_f, 1, memoized_f, ∇f) + + @constraint(upper_model, y_star == op_f(x)) + + @time optimize!(upper_model) + + @test objective_value(upper_model) ≈ 15.85 atol=0.05 + @test value(x) ≈ 2.0 atol=0.05 + @test value(y) ≈ 3.0 atol=0.05 +end + +################################################ +#= +# Test Stratigic Bidding +=# +################################################ + +function build_lower() + lower_model = Model(Ipopt.Optimizer) + + @variable(lower_model, qS_p ∈ MOI.Parameter(1.0)) + @variable(lower_model, 0 <= gS <= 100) + @variable(lower_model, 0 <= gR1 <= 40) + @variable(lower_model, 0 <= gR2 <= 40) + @variable(lower_model, 0 <= gD <= 100) + @objective(lower_model, Min, 0.5gR1^2 + 1gR2^2 + 2gD^2) + @constraint(lower_model, bid, gS <= qS_p) + @constraint(lower_model, demand_equilibrium, gS + gR1 + gR2 + gD == 100) + primal_vars = [gS; gR1; gR2; gD] + params = [qS_p] + return lower_model, qS_p, demand_equilibrium, primal_vars, params, bid +end + +function build_upper() + upper_model = Model(Ipopt.Optimizer) + + @variable(upper_model, 0 <= qS <= 100) + @variable(upper_model, lambda) + @objective(upper_model, Max, lambda*qS) + return upper_model, lambda, qS +end + +function test_bilevel_strategic_bidding() + # test derivative of the dual of the demand equilibrium constraint + lower_model, qS_p, demand_equilibrium, primal_vars, params, bid = build_lower() + upper_model, lambda, qS = build_upper() + + optimize!(lower_model) + + evaluator, cons = create_evaluator(lower_model; x=[primal_vars; params]) + Δs, sp = compute_sensitivity(evaluator, cons, [0.5]; primal_vars=primal_vars, params=params) + + set_parameter_value(qS_p, 1.5) + + optimize!(lower_model) + + @test dual(demand_equilibrium) ≈ sp[6] + @test dual(bid) ≈ sp[7] + + # test bilevel strategic bidding + + lower_model, qS_p, demand_equilibrium, primal_vars, params = build_lower() + upper_model, lambda, qS = build_upper() + + evaluator, cons = create_evaluator(lower_model; x=[primal_vars; params]) + + function f(qS_val) + set_parameter_value(qS_p, qS_val) + optimize!(lower_model) + @assert is_solved_and_feasible(lower_model) + return dual(demand_equilibrium) + end + + function ∇f(qS_val) + @assert value(qS_p) == qS_val + Δs, sp = compute_sensitivity(evaluator, cons, [0.001]; primal_vars=primal_vars, params=params) + return Δs[6] + end + + memoized_f = memoize(f) + + @operator(upper_model, op_f, 1, memoized_f, ∇f) + + @constraint(upper_model, lambda == op_f(qS)) + + # 0.68s + optimize!(upper_model) + @test solve_time(upper_model) < 1.0 + @test objective_value(upper_model) ≈ 0.00508 rtol=1e-2 + + ###### no derivative + + lower_model, qS_p, demand_equilibrium, primal_vars, params = build_lower() + upper_model, lambda, qS = build_upper() + + evaluator, cons = create_evaluator(lower_model; x=[primal_vars; params]) + + function f(qS_val) + set_parameter_value(qS_p, qS_val) + optimize!(lower_model) + @assert is_solved_and_feasible(lower_model) + return dual(demand_equilibrium) + end + + # function ∇f(qS_val) + # @assert value(qS_p) == qS_val + # Δs, sp = compute_sensitivity(evaluator, cons, [0.00]; primal_vars=primal_vars, params=params) + # return Δs[6] + # end + + memoized_f = memoize(f) + + @operator(upper_model, op_f, 1, memoized_f, (x) -> 0.0) + + @constraint(upper_model, lambda == op_f(qS)) + + # 1.09s + optimize!(upper_model) + @test solve_time(upper_model) > 1.0 + @test objective_value(upper_model) ≈ 0.00508 rtol=1e-2 +end diff --git a/moonlanding.jl b/moonlanding.jl new file mode 100644 index 00000000..27ad2771 --- /dev/null +++ b/moonlanding.jl @@ -0,0 +1,75 @@ +function moonlander_JMP(;target::Array{Float64}=[5.0, 5.0],nh::Int64=500, _I = 0.1) + ## parameters + if size(target) != (2,) + error("The input matrix must be 2x1.") + end + m = 1.0 + g = 9.81 + D = 1.0 + max_thrust = 2*g + + ## define the problem + model = JuMP.Model() + + @variable(model, Isp ∈ MOI.Parameter(_I)) + + @variables(model, begin + 0.0 <= tf + # state variables + p1[k=0:nh] + p2[k=0:nh] + dp1[k=0:nh] + dp2[k=0:nh] + theta[k=0:nh] + dtheta[k=0:nh] + + # control variables + 0 <= F1[k=0:nh] <= max_thrust, (start = 5.0) + 0 <= F2[k=0:nh] <= max_thrust, (start = 5.0) + end) + + # Initial and final conditions + @constraints(model, begin + p1[0] == 0.0 + p2[0] == 0.0 + dp1[0] == 0.0 + dp2[0] == 0.0 + theta[0] == 0.0 + dtheta[0] == 0.0 + p1[nh] == target[1] + p2[nh] == target[2] + dp1[nh] == 0.0 + dp2[nh] == 0.0 + end) + + #dynamics + @expressions(model, begin + F_r[k=0:nh], [cos(theta[k]) -sin(theta[k]) p1[k]; + sin(theta[k]) cos(theta[k]) p2[k]; + 0.0 0.0 1.0] + end) + @expressions(model, begin + F_tot[k=0:nh], (F_r[k] * [0; F1[k] + F2[k]; 0])[1:2] + end) + @expressions(model, begin + ddp1[k=0:nh], (1/m) * F_tot[k][1] + ddp2[k=0:nh], (1/m) * F_tot[k][2] - g + ddtheta[k=0:nh], (1/Isp) * (D/2) * (F2[k] - F1[k]) + end) + @expressions(model, begin + step, tf / nh + end) + + @constraints(model, begin + d_p1[k=1:nh], p1[k] == p1[k-1] + 0.5 * step * (dp1[k] + dp1[k-1]) + d_p2[k=1:nh], p2[k] == p2[k-1] + 0.5 * step * (dp2[k] + dp2[k-1]) + d_dp1[k=1:nh], dp1[k] == dp1[k-1] + 0.5 * step * (ddp1[k] + ddp1[k-1]) + d_dp2[k=1:nh], dp2[k] == dp2[k-1] + 0.5 * step * (ddp2[k] + ddp2[k-1]) + d_theta[k=1:nh], theta[k] == theta[k-1] + 0.5 * step * (dtheta[k] + dtheta[k-1]) + d_dtheta[k=1:nh], dtheta[k] == dtheta[k-1] + 0.5 * step * (ddtheta[k] + ddtheta[k-1]) + end) + + @objective(model, Min, tf) + + return model, tf, Isp +end \ No newline at end of file diff --git a/nlp_utilities.jl b/nlp_utilities.jl new file mode 100644 index 00000000..78d2c0d4 --- /dev/null +++ b/nlp_utilities.jl @@ -0,0 +1,596 @@ +using JuMP +using MathOptInterface +import MathOptInterface: ConstraintSet, CanonicalConstraintFunction +using SparseArrays +using LinearAlgebra +import JuMP.index + +""" + create_nlp_model(model::JuMP.Model) + +Create a Nonlinear Programming (NLP) model from a JuMP model. +""" +function create_nlp_model(model::JuMP.Model) + rows = Vector{ConstraintRef}(undef, 0) + nlp = MOI.Nonlinear.Model() + for (F, S) in list_of_constraint_types(model) + if F <: VariableRef && !(S <: MathOptInterface.EqualTo{Float64}) + continue # Skip variable bounds + end + for ci in all_constraints(model, F, S) + push!(rows, ci) + object = constraint_object(ci) + MOI.Nonlinear.add_constraint(nlp, object.func, object.set) + end + end + MOI.Nonlinear.set_objective(nlp, objective_function(model)) + return nlp, rows +end + +""" + fill_off_diagonal(H) + +Filling the off-diagonal elements of a sparse matrix to make it symmetric. +""" +function fill_off_diagonal(H) + ret = H + H' + row_vals = SparseArrays.rowvals(ret) + non_zeros = SparseArrays.nonzeros(ret) + for col in 1:size(ret, 2) + for i in SparseArrays.nzrange(ret, col) + if col == row_vals[i] + non_zeros[i] /= 2 + end + end + end + return ret +end + +sense_mult(x) = JuMP.objective_sense(owner_model(x)) == MOI.MIN_SENSE ? 1.0 : -1.0 +sense_mult(x::Vector) = sense_mult(x[1]) + +""" + compute_optimal_hessian(evaluator::MOI.Nonlinear.Evaluator, rows::Vector{ConstraintRef}, x::Vector{VariableRef}) + +Compute the optimal Hessian of the Lagrangian. +""" +function compute_optimal_hessian(evaluator::MOI.Nonlinear.Evaluator, rows::Vector{ConstraintRef}, x::Vector{VariableRef}) + sense_multiplier = objective_sense(owner_model(x[1])) == MOI.MIN_SENSE ? 1.0 : -1.0 + hessian_sparsity = MOI.hessian_lagrangian_structure(evaluator) + I = [i for (i, _) in hessian_sparsity] + J = [j for (_, j) in hessian_sparsity] + V = zeros(length(hessian_sparsity)) + # The signals are being sdjusted to match the Ipopt convention (inner.mult_g) + # but we don't know if we need to adjust the objective function multiplier + MOI.eval_hessian_lagrangian(evaluator, V, value.(x), 1.0, - sense_multiplier * dual.(rows)) + H = SparseArrays.sparse(I, J, V, length(x), length(x)) + return fill_off_diagonal(H) +end + +""" + compute_optimal_jacobian(evaluator::MOI.Nonlinear.Evaluator, rows::Vector{ConstraintRef}, x::Vector{VariableRef}) + +Compute the optimal Jacobian of the constraints. +""" +function compute_optimal_jacobian(evaluator::MOI.Nonlinear.Evaluator, rows::Vector{ConstraintRef}, x::Vector{VariableRef}) + jacobian_sparsity = MOI.jacobian_structure(evaluator) + I = [i for (i, _) in jacobian_sparsity] + J = [j for (_, j) in jacobian_sparsity] + V = zeros(length(jacobian_sparsity)) + MOI.eval_constraint_jacobian(evaluator, V, value.(x)) + A = SparseArrays.sparse(I, J, V, length(rows), length(x)) + return A +end + +""" + compute_optimal_hess_jac(evaluator::MOI.Nonlinear.Evaluator, rows::Vector{ConstraintRef}, x::Vector{VariableRef}) + +Compute the optimal Hessian of the Lagrangian and Jacobian of the constraints. +""" +function compute_optimal_hess_jac(evaluator::MOI.Nonlinear.Evaluator, rows::Vector{ConstraintRef}, x::Vector{VariableRef}) + hessian = compute_optimal_hessian(evaluator, rows, x) + jacobian = compute_optimal_jacobian(evaluator, rows, x) + + return hessian, jacobian +end + +""" + all_primal_vars(model::Model) + +Get all the primal variables in the model. +""" +all_primal_vars(model::Model) = filter(x -> !is_parameter(x), all_variables(model)) + +""" + all_params(model::Model) + +Get all the parameters in the model. +""" +all_params(model::Model) = filter(x -> is_parameter(x), all_variables(model)) + +""" + create_evaluator(model::Model; x=all_variables(model)) + +Create an evaluator for the model. +""" +JuMP.index(x::JuMP.Containers.DenseAxisArray) = index.(x).data + +function create_evaluator(model::Model; x=all_variables(model)) + nlp, rows = create_nlp_model(model) + backend = MOI.Nonlinear.SparseReverseMode() + evaluator = MOI.Nonlinear.Evaluator(nlp, backend, vcat(index.(x)...)) + MOI.initialize(evaluator, [:Hess, :Jac]) + return evaluator, rows +end + +""" + is_inequality(con::ConstraintRef) + +Check if the constraint is an inequality. +""" +# function is_inequality(con::ConstraintRef) +# set_type = typeof(MOI.get(owner_model(con), MOI.ConstraintSet(), con)) +# return set_type <: MOI.LessThan || set_type <: MOI.GreaterThan +# end + +function is_less_inequality(con::ConstraintRef) + set_type = typeof(MOI.get(owner_model(con), MOI.ConstraintSet(), con)) + return set_type <: MOI.LessThan +end + +function is_greater_inequality(con::ConstraintRef) + set_type = typeof(MOI.get(owner_model(con), MOI.ConstraintSet(), con)) + return set_type <: MOI.GreaterThan +end + +function is_inequality(con::ConstraintRef) + return is_less_inequality(con) || is_greater_inequality(con) +end + +""" + find_inequealities(cons::Vector{ConstraintRef}) + +Find the indices of the inequality constraints. +""" +function find_inequealities(cons::Vector{C}) where C<:ConstraintRef + leq_locations = zeros(length(cons)) + geq_locations = zeros(length(cons)) + for i in 1:length(cons) + if is_less_inequality(cons[i]) + leq_locations[i] = true + end + if is_greater_inequality(cons[i]) + geq_locations[i] = true + end + end + return findall(x -> x ==1, leq_locations), findall(x -> x ==1, geq_locations) +end + +""" + get_slack_inequality(con::ConstraintRef) + +Get the reference to the canonical function that is equivalent to the slack variable of the inequality constraint. +""" +function get_slack_inequality(con::ConstraintRef) + set_type = typeof(MOI.get(owner_model(con), MOI.ConstraintSet(), con)) + obj = constraint_object(con) + if set_type <: MOI.LessThan + # c(x) <= b --> slack = c(x) - b | slack <= 0 + return obj.func - obj.set.upper + end + # c(x) >= b --> slack = c(x) - b | slack >= 0 + return obj.func - obj.set.lower +end + +""" + compute_solution_and_bounds(primal_vars::Vector{VariableRef}, cons::Vector{C}) where C<:ConstraintRef + +Compute the solution and bounds of the primal variables. +""" +function compute_solution_and_bounds(primal_vars::Vector{VariableRef}, cons::Vector{C}) where {C<:ConstraintRef} + sense_multiplier = sense_mult(primal_vars) + num_vars = length(primal_vars) + leq_locations, geq_locations = find_inequealities(cons) + ineq_locations = vcat(geq_locations, leq_locations) + num_leq = length(leq_locations) + num_geq = length(geq_locations) + num_ineq = num_leq + num_geq + slack_vars = [get_slack_inequality(cons[i]) for i in ineq_locations] + has_up = findall(x -> has_upper_bound(x), primal_vars) + has_low = findall(x -> has_lower_bound(x), primal_vars) + + # Primal solution + X = value.([primal_vars; slack_vars]) + + # value and dual of the lower bounds + V_L = spzeros(num_vars+num_ineq) + X_L = spzeros(num_vars+num_ineq) + for (i, j) in enumerate(has_low) + V_L[i] = dual.(LowerBoundRef(primal_vars[j])) * sense_multiplier + # + if sense_multiplier == 1.0 + V_L[i] <= -1e-6 && @info "Dual of lower bound must be positive" i V_L[i] + else + V_L[i] >= 1e-6 && @info "Dual of lower bound must be negative" i V_L[i] + end + # + X_L[i] = JuMP.lower_bound(primal_vars[j]) + end + for (i, con) in enumerate(cons[geq_locations]) + # By convention jump dual will allways be positive for geq constraints + # but for ipopt it will be positive if min problem and negative if max problem + V_L[num_vars+i] = dual.(con) * (sense_multiplier) + # + if sense_multiplier == 1.0 + V_L[num_vars+i] <= -1e-6 && @info "Dual of geq constraint must be positive" i V_L[num_vars+i] + else + V_L[num_vars+i] >= 1e-6 && @info "Dual of geq constraint must be negative" i V_L[num_vars+i] + end + end + # value and dual of the upper bounds + V_U = spzeros(num_vars+num_ineq) + X_U = spzeros(num_vars+num_ineq) + for (i, j) in enumerate(has_up) + V_U[i] = dual.(UpperBoundRef(primal_vars[j])) * (- sense_multiplier) + # + if sense_multiplier == 1.0 + V_U[i] <= -1e-6 && @info "Dual of upper bound must be positive" i V_U[i] + else + V_U[i] >= 1e-6 && @info "Dual of upper bound must be negative" i V_U[i] + end + # + X_U[i] = JuMP.upper_bound(primal_vars[j]) + end + for (i, con) in enumerate(cons[leq_locations]) + # By convention jump dual will allways be negative for leq constraints + # but for ipopt it will be positive if min problem and negative if max problem + V_U[num_vars+i] = dual.(con) * (- sense_multiplier) + # + if sense_multiplier == 1.0 + V_U[num_vars+i] <= -1e-6 && @info "Dual of leq constraint must be positive" i V_U[num_vars+i] + else + V_U[num_vars+i] >= 1e-6 && @info "Dual of leq constraint must be negative" i V_U[num_vars+i] + end + end + return X, V_L, X_L, V_U, X_U, leq_locations, geq_locations, ineq_locations, vcat(has_up, collect(num_vars+num_geq+1:num_vars+num_geq+num_leq)), vcat(has_low, collect(num_vars+1:num_vars+num_geq)) +end + +""" + build_M_N(evaluator::MOI.Nonlinear.Evaluator, cons::Vector{ConstraintRef}, + primal_vars::Vector{VariableRef}, params::Vector{VariableRef}, + _X::Vector, _V_L::Vector, _X_L::Vector, _V_U::Vector, _X_U::Vector, ineq_locations::Vector{Z}, + has_up::Vector{Z}, has_low::Vector{Z} +) where {Z<:Integer} + +Build the M (KKT Jacobian w.r.t. solution) and N (KKT Jacobian w.r.t. parameters) matrices for the sensitivity analysis. +""" +function build_M_N(evaluator::MOI.Nonlinear.Evaluator, cons::Vector{ConstraintRef}, + primal_vars::Vector{VariableRef}, params::Vector{VariableRef}, + _X::AbstractVector, _V_L::AbstractVector, _X_L::AbstractVector, _V_U::AbstractVector, _X_U::AbstractVector, leq_locations::Vector{Z}, geq_locations::Vector{Z}, ineq_locations::Vector{Z}, + has_up::Vector{Z}, has_low::Vector{Z} +) where {Z<:Integer} + @assert all(x -> is_parameter(x), params) "All parameters must be parameters" + + # Setting + num_vars = length(primal_vars) + num_parms = length(params) + num_cons = length(cons) + num_ineq = length(ineq_locations) + all_vars = [primal_vars; params] + num_low = length(has_low) + num_up = length(has_up) + + # Primal solution + X_lb = spzeros(num_low, num_low) + X_ub = spzeros(num_up, num_up) + V_L = spzeros(num_low, num_vars + num_ineq) + V_U = spzeros(num_up, num_vars + num_ineq) + I_L = spzeros(num_vars + num_ineq, num_low) + I_U = spzeros(num_vars + num_ineq, num_up) + + # value and dual of the lower bounds + for (i, j) in enumerate(has_low) + V_L[i, j] = _V_L[j] + X_lb[i, i] = _X[j] - _X_L[j] + I_L[j, i] = -1 + end + # value and dual of the upper bounds + for (i, j) in enumerate(has_up) + V_U[i, j] = _V_U[j] + X_ub[i, i] = _X_U[j] - _X[j] + I_U[j, i] = 1 + end + + # Function Derivatives + hessian, jacobian = compute_optimal_hess_jac(evaluator, cons, all_vars) + + # Hessian of the lagrangian wrt the primal variables + W = spzeros(num_vars + num_ineq, num_vars + num_ineq) + W[1:num_vars, 1:num_vars] = hessian[1:num_vars, 1:num_vars] + # Jacobian of the constraints + A = spzeros(num_cons, num_vars + num_ineq) + # A is the Jacobian of: c(x) = b and c(x) <= b and c(x) >= b, possibly all mixed up. + # Each of the will be re-written as: + # c(x) - b = 0 + # c(x) - b - su = 0, su <= 0 + # c(x) - b - sl = 0, sl >= 0 + # Jacobian of the constraints wrt the primal variables + A[:, 1:num_vars] = jacobian[:, 1:num_vars] + # Jacobian of the constraints wrt the slack variables + for (i,j) in enumerate(geq_locations) + A[j, num_vars+i] = -1 + end + for (i,j) in enumerate(leq_locations) + A[j, num_vars+i] = -1 + end + # Partial second derivative of the lagrangian wrt primal solution and parameters + ∇ₓₚL = spzeros(num_vars + num_ineq, num_parms) + ∇ₓₚL[1:num_vars, :] = hessian[1:num_vars, num_vars+1:end] + # Partial derivative of the equality constraintswith wrt parameters + ∇ₚC = jacobian[:, num_vars+1:end] + + # M matrix + # M = [ + # [W A' -I I]; + # [A 0 0 0]; + # [V_L 0 (X - X_L) 0] + # [V_U 0 0 0 (X_U - X)] + # ] + len_w = num_vars + num_ineq + M = spzeros(len_w + num_cons + num_low + num_up, len_w + num_cons + num_low + num_up) + + M[1:len_w, 1:len_w] = W + M[1:len_w, len_w + 1 : len_w + num_cons] = A' + M[len_w+1:len_w+num_cons, 1:len_w] = A + M[1:len_w, len_w+num_cons+1:len_w+num_cons+num_low] = I_L + M[len_w+num_cons+1:len_w+num_cons+num_low, 1:len_w] = V_L + M[len_w+num_cons+1:len_w+num_cons+num_low, len_w+num_cons+1:len_w+num_cons+num_low] = X_lb + M[len_w+num_cons+num_low+1:len_w+num_cons+num_low+num_up, 1:len_w] = V_U + M[len_w+num_cons+num_low+1:len_w+num_cons+num_low+num_up, len_w+num_cons+num_low+1:len_w+num_cons+num_low+num_up] = X_ub + M[1:len_w, len_w+num_cons+num_low+1:end] = I_U + + # N matrix + # N = [∇ₓₚL ; ∇ₚC; zeros(num_low + num_up, num_parms)] + N = spzeros(len_w + num_cons + num_low + num_up, num_parms) + N[1:len_w, :] = ∇ₓₚL + N[len_w+1:len_w+num_cons, :] = ∇ₚC + + return M, N +end + +function inertia_corrector_factorization(M::SparseMatrixCSC, num_w, num_cons; st=1e-6, max_corrections=50) + # Factorization + K = lu(M; check=false) + # Inertia correction + status = K.status + num_c = 0 + diag_mat = ones(size(M, 1)) + diag_mat[num_w+1:num_w+num_cons] .= -1 + diag_mat = sparse(diagm(diag_mat)) + while status == 1 && num_c < max_corrections + println("Inertia correction") + M = M + st * diag_mat + K = lu(M; check=false) + status = K.status + num_c += 1 + end + if status != 0 + @warn "Inertia correction failed" + return nothing + end + return K +end + +function inertia_corrector_factorization(M; st=1e-6, max_corrections=50) + num_c = 0 + if cond(M) > 1/st + @warn "Inertia correction" + M = M + st * I(size(M, 1)) + num_c += 1 + end + while cond(M) > 1/st && num_c < max_corrections + M = M + st * I(size(M, 1)) + num_c += 1 + end + if num_c == max_corrections + @warn "Inertia correction failed" + return nothing + end + return lu(M) +end + +""" + compute_derivatives_no_relax(evaluator::MOI.Nonlinear.Evaluator, cons::Vector{ConstraintRef}, + primal_vars::Vector{VariableRef}, params::Vector{VariableRef}, + _X::Vector, _V_L::Vector, _X_L::Vector, _V_U::Vector, _X_U::Vector, ineq_locations::Vector{Z}, + has_up::Vector{Z}, has_low::Vector{Z} + ) where {Z<:Integer} + +Compute the derivatives of the solution w.r.t. the parameters without accounting for active set changes. +""" +function compute_derivatives_no_relax(evaluator::MOI.Nonlinear.Evaluator, cons::Vector{ConstraintRef}, + primal_vars::Vector{VariableRef}, params::Vector{VariableRef}, + _X::AbstractVector, _V_L::AbstractVector, _X_L::AbstractVector, _V_U::AbstractVector, _X_U::AbstractVector, leq_locations::Vector{Z}, geq_locations::Vector{Z}, ineq_locations::Vector{Z}, + has_up::Vector{Z}, has_low::Vector{Z} +) where {Z<:Integer} + num_bounds = length(has_up) + length(has_low) + M, N = build_M_N(evaluator, cons, primal_vars, params, _X, _V_L, _X_L, _V_U, _X_U, leq_locations, geq_locations, ineq_locations, has_up, has_low) + + # Sesitivity of the solution (primal-dual_constraints-dual_bounds) w.r.t. the parameters + K = inertia_corrector_factorization(M, length(primal_vars) + length(ineq_locations), length(cons)) # Factorization + if isnothing(K) + return zeros(size(M, 1), size(N, 2)), K, N + end + ∂s = zeros(size(M, 1), size(N, 2)) + # ∂s = - (K \ N) # Sensitivity + ldiv!(∂s, K, N) + ∂s = - ∂s + + return ∂s, K, N +end + +""" + fix_and_relax(E, K, N, r1, Δp, num_w, num_cons) + +Fix the violations and relax complementary slackness. +""" +function fix_and_relax(E, K, N, r1, Δp) + rs = (N * Δp)[:,:] + # C = −E' inv(K) E + # C = - E' * (K \ E) + C = zeros(size(K, 1), size(E, 2)) + ldiv!(C, K, E) + C = - E' * C + # C ∆ν¯ = E' inv(K) rs − r1 + # ∆ν¯ = C \ (E' * (K \ rs) - r1) + aux = zeros(size(K, 1), size(rs, 2)) + ldiv!(aux, K, rs) + C_fact = inertia_corrector_factorization(C) + if isnothing(C_fact) + return zeros(size(K, 1), size(aux, 2)) + end + ∆ν¯ = C_fact \ (E' * aux - r1) + # K ∆s = − (rs + E∆ν¯) + # ∆s = K \ (- (rs + E * ∆ν¯)) + aux = - (rs + E * ∆ν¯)[:,:] + ∆s = zeros(size(K, 1), size(aux, 2)) + ldiv!(∆s, K, aux) + + return ∆s +end + +""" + approximate_solution(X, Λ, V_L, V_U, Δs) + +First order approximation the primal-dual solution. +""" +function approximate_solution(X::Vector{T}, Λ, V_L, V_U, Δs) where T + num_var = length(X) + num_cons = length(Λ) + num_low = length(V_L) + num_up = length(V_U) + sp = Vector{T}(undef, num_var + num_cons + num_low + num_up) + sp[1:num_var] = X + sp[num_var+1:num_var+num_cons] = Λ + sp[num_var+num_cons+1:num_var+num_cons+num_low] = V_L + sp[num_var+num_cons+num_low+1:num_var+num_cons+num_low+num_up] = V_U + sp = sp + Δs + + return sp +end + +""" + find_violations(X, sp, X_L, X_U, V_U, V_L, has_up, has_low, num_cons, tol) + +Find the bound violations of the primal-dual solution first order approximation based on the (unrelaxed) sensitivity estimation. +""" +function find_violations(X, sp, X_L, X_U, V_U, V_L, has_up, has_low, num_cons, tol, ismin) + num_w = length(X) + num_low = length(has_low) + num_up = length(has_up) + + _E = [] + r1 = [] + for (j, i) in enumerate(has_low) + if sp[i] < X_L[i] - tol + println("Violation LB Var: ", i, " ", sp[i], " ", X_L[i]) + push!(_E, i) + push!(r1, X[i] - X_L[i]) + end + if ismin + if sp[num_w+num_cons+j] < -tol + println("Violation LB Dual: ", i, " ", sp[num_w+num_cons+j], " ", V_L[i]) + push!(_E, num_w+num_cons+j) + push!(r1, V_L[i]) + end + else + if sp[num_w+num_cons+j] > tol + println("Violation LB Dual: ", i, " ", sp[num_w+num_cons+j], " ", V_L[i]) + push!(_E, num_w+num_cons+j) + push!(r1, V_L[i]) + end + end + end + for (j, i) in enumerate(has_up) + if sp[i] > X_U[i] + tol + println("Violation UB Var: ", i, " ", sp[i], " ", X_U[i]) + push!(_E, i) + push!(r1, X[i] - X_U[i]) + end + if ismin + if sp[num_w+num_cons+num_low+j] < -tol + println("Violation UB Dual: ", i, " ", sp[num_w+num_cons+num_low+j], " ", V_U[i]) + push!(_E, num_w+num_cons+num_low+j) + push!(r1, V_U[i]) + end + else + if sp[num_w+num_cons+num_low+j] > tol + println("Violation UB Dual: ", i, " ", sp[num_w+num_cons+num_low+j], " ", V_U[i]) + push!(_E, num_w+num_cons+num_low+j) + push!(r1, V_U[i]) + end + end + end + + E = spzeros(num_w + num_cons + num_low + num_up, length(_E)) + for (i, j) in enumerate(_E) + E[j, i] = 1 + end + + return E, r1 +end + + +function compute_sensitivity(evaluator::MOI.Nonlinear.Evaluator, cons::Vector{ConstraintRef}, + Δp; primal_vars=all_primal_vars(model), params=all_params(model), tol=1e-6 +) + ismin = sense_mult(primal_vars) == 1.0 + sense_multiplier = sense_mult(primal_vars) + num_cons = length(cons) + num_var = length(primal_vars) + # Solution and bounds + X, V_L, X_L, V_U, X_U, leq_locations, geq_locations, ineq_locations, has_up, has_low = compute_solution_and_bounds(primal_vars, cons) + # Compute derivatives + num_w = length(X) + num_lower = length(has_low) + # ∂s = [∂x; ∂λ; ∂ν_L; ∂ν_U] + ∂s, K, N = compute_derivatives_no_relax(evaluator, cons, primal_vars, params, X, V_L, X_L, V_U, X_U, leq_locations, geq_locations, ineq_locations, has_up, has_low) + if isnothing(Δp) + ## Adjust signs based on JuMP convention + # Duals + ∂s[num_w+1:num_w+num_cons, :] *= -sense_multiplier + # Dual bounds lower + ∂s[num_w+num_cons+1:num_w+num_cons+num_lower, :] *= sense_multiplier + # Dual bounds upper + ∂s[num_w+num_cons+num_lower+1:end, :] *= -sense_multiplier + sp = approximate_solution(X, dual.(cons), V_L[has_low] .* (sense_multiplier), V_U[has_up].* (-sense_multiplier), ∂s * ones(size(∂s, 2))) + return ∂s, sp + end + Δs = ∂s * Δp + sp = approximate_solution(X, - sense_multiplier * dual.(cons), V_L[has_low], V_U[has_up], Δs) + # Linearly appoximated solution + E, r1 = find_violations(X, sp, X_L, X_U, V_U, V_L, has_up, has_low, num_cons, tol, ismin) + if !isempty(r1) + @warn "Relaxation needed" + Δs = fix_and_relax(E, K, N, r1, Δp) + end + ## Adjust signs based on JuMP convention + # Duals + Δs[num_w+1:num_w+num_cons, :] *= -sense_multiplier + # Dual bounds lower + Δs[num_w+num_cons+1:num_w+num_cons+num_lower, :] *= sense_multiplier + # Dual bounds upper + Δs[num_w+num_cons+num_lower+1:end, :] *= -sense_multiplier + sp = approximate_solution(X, dual.(cons), V_L[has_low] .* (sense_multiplier), V_U[has_up].* (-sense_multiplier), Δs) + return Δs, sp +end + +""" + compute_sensitivity(model::Model; primal_vars=all_primal_vars(model), params=all_params(model)) + +Compute the sensitivity of the solution given sensitivity of the parameters (Δp). +""" +function compute_sensitivity(model::Model, Δp; primal_vars=all_primal_vars(model), params=all_params(model)) + evaluator, cons = create_evaluator(model; x=[primal_vars; params]) + return compute_sensitivity(evaluator, cons, Δp; primal_vars=primal_vars, params=params), evaluator, cons +end diff --git a/nlp_utilities_test.jl b/nlp_utilities_test.jl new file mode 100644 index 00000000..a76426d5 --- /dev/null +++ b/nlp_utilities_test.jl @@ -0,0 +1,299 @@ +using JuMP +using Ipopt +using Test +using FiniteDiff + +include("nlp_utilities.jl") + +include("problem_list.jl") + +################################################ +#= +# Test JuMP Hessian and Jacobian + +From JuMP Tutorial for Querying Hessians: +https://github.com/jump-dev/JuMP.jl/blob/301d46e81cb66c74c6e22cd89fb89ced740f157b/docs/src/tutorials/nonlinear/querying_hessians.jl#L67-L72 +=# +################################################ + +function analytic_hessian(x, σ, μ, p) + g_1_H = [2.0 0.0; 0.0 0.0] + g_2_H = p[1] * [2.0 2.0; 2.0 2.0] + f_H = zeros(2, 2) + f_H[1, 1] = 2.0 + p[3] * 12.0 * x[1]^2 - p[3] * 4.0 * x[2] + f_H[1, 2] = f_H[2, 1] = -p[3] * 4.0 * x[1] + f_H[2, 2] = p[3] * 2.0 + return σ * f_H + μ' * [g_1_H, g_2_H] +end + +function analytic_jacobian(x, p) + g_1_J = [ + 2.0 * x[1], # ∂g_1/∂x_1 + 0.0, # ∂g_1/∂x_2 + -1.0, # ∂g_1/∂p_1 + 0.0, # ∂g_1/∂p_2 + 0.0 # ∂g_1/∂p_3 + ] + g_2_J = [ + p[1] * 2.0 * (x[1] + x[2]), # ∂g_2/∂x_1 + 2.0 * (x[1] + x[2]), # ∂g_2/∂x_2 + (x[1] + x[2])^2, # ∂g_2/∂p_1 + -1.0, # ∂g_2/∂p_2 + 0.0 # ∂g_2/∂p_3 + ] + return hcat(g_2_J, g_1_J)'[:,:] +end + +function test_create_evaluator(model, x) + @testset "Create NLP model" begin + nlp, rows = create_nlp_model(model) + @test nlp isa MOI.Nonlinear.Model + @test rows isa Vector{ConstraintRef} + end + @testset "Create Evaluator" begin + evaluator, rows = create_evaluator(model; x = x) + @test evaluator isa MOI.Nonlinear.Evaluator + @test rows isa Vector{ConstraintRef} + end +end + +function test_compute_optimal_hess_jacobian() + @testset "Compute Optimal Hessian and Jacobian" begin + # Model + model, x, cons, params = create_nonlinear_jump_model() + # Optimize + optimize!(model) + @assert is_solved_and_feasible(model) + # Create evaluator + test_create_evaluator(model, [x; params]) + evaluator, rows = create_evaluator(model; x = [x; params]) + # Compute Hessian and Jacobian + num_var = length(x) + full_hessian, full_jacobian = compute_optimal_hess_jac(evaluator, rows, [x; params]) + hessian = full_hessian[1:num_var, 1:num_var] + # Check Hessian + @test all(hessian .≈ analytic_hessian(value.(x), 1.0, -dual.(cons), value.(params))) + # TODO: Test hessial of parameters + # Check Jacobian + @test all(full_jacobian .≈ analytic_jacobian(value.(x), value.(params))) + end +end + +################################################ +#= +# Test Derivatives and Sensitivity: QP problem + +From sIpopt paper: https://optimization-online.org/2011/04/3008/ +=# +################################################ + +function test_compute_derivatives() + @testset "Compute Derivatives No Inequalities: ismin=$ismin" for ismin in [true, false] + # Model + # sign_fix = ismin ? 1.0 : -1.0 + model, primal_vars, cons, params = create_nonlinear_jump_model_sipopt(;ismin=ismin) + optimize!(model) + @assert is_solved_and_feasible(model) + # Analytical solutions case b + pb = [4.5, 1.0] + s_pb = [0.5, 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0] + @assert all(isapprox.([value.(primal_vars); dual.(cons); dual.(LowerBoundRef.(primal_vars))], s_pb; atol = 1e-6)) + # Analytical solutions case a + pa = [5.0, 1.0] + s_pa = [0.6327, 0.3878, 0.0204, 0.1633, 0.2857, 0, 0, 0] + set_parameter_value.(params, pa) + optimize!(model) + @assert is_solved_and_feasible(model) + @assert all(isapprox.([value.(primal_vars); dual.(cons); dual.(LowerBoundRef.(primal_vars))], s_pa; atol = 1e-4)) + # Compute derivatives without accounting for active set changes + evaluator, rows = create_evaluator(model; x=[primal_vars; params]) + ∂s, s_pb_approx = compute_sensitivity(evaluator, rows, nothing; primal_vars, params) + # Check linear approx s_pb + Δp = pb - pa + s_pb_approx_violated = s_pa + ∂s * Δp + @test all(isapprox.([0.5765; 0.3775; -0.0459; 0.1327; 0.3571; 0.0; 0.0; 0.0], s_pb_approx_violated; atol = 1e-4)) + # Account for active set changes + Δs, s_pb_approx = compute_sensitivity(evaluator, rows, Δp; primal_vars, params) + # s_pb_approx = s_pa .+ Δs + @test all(isapprox.([0.5, 0.5, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0], s_pb_approx; atol = 1e-6)) + end +end + +################################################ +#= +# Test Sensitivity through analytical +=# +################################################ + + +# f(x, p) = 0 +# x = g(p) +# ∂x/∂p = ∂g/∂p + +DICT_PROBLEMS_Analytical_no_cc = Dict( + "geq no impact" => (p_a=[1.5], Δp=[0.2], Δs_a=[0.0; -0.2; 0.0; 0.0; 0.0; 0.0; 0.0], model_generator=create_jump_model_1), + "geq impact" => (p_a=[2.1], Δp=[0.2], Δs_a=[0.2; 0.0; 0.2; 0.4; 0.0; 0.4; 0.0], model_generator=create_jump_model_1), + "geq bound impact" => (p_a=[2.1], Δp=[0.2], Δs_a=[0.2; 0.0; 0.4; 0.0; 0.4], model_generator=create_jump_model_2), + "leq no impact" => (p_a=[-1.5], Δp=[-0.2], Δs_a=[0.0; 0.2; 0.0; 0.0; 0.0; 0.0; 0.0], model_generator=create_jump_model_3), + "leq impact" => (p_a=[-2.1], Δp=[-0.2], Δs_a=[-0.2; 0.0; -0.2], model_generator=create_jump_model_3), + "leq no impact max" => (p_a=[2.1], Δp=[0.2], Δs_a=[0.0; -0.2; 0.0; 0.0; 0.0], model_generator=create_jump_model_4), + "leq impact max" => (p_a=[1.5], Δp=[0.2], Δs_a=[0.2; 0.0; 0.2], model_generator=create_jump_model_4), + "geq no impact max" => (p_a=[1.5], Δp=[0.2], Δs_a=[0.0; -0.2; 0.0; 0.0; 0.0], model_generator=create_jump_model_5), + "geq impact max" => (p_a=[2.1], Δp=[0.2], Δs_a=[0.2; 0.0; 0.2], model_generator=create_jump_model_5), + # "softmax" => (p_a=collect(1.0:0.1:2.0), Δp=collect(-0.1:0.02:0.1), Δs_a=jacobian(softmax, p_a)[1] * Δp, model_generator=create_jump_model_6) +) + +DICT_PROBLEMS_Analytical_cc = Dict( + "geq active constraint change" => (p_a=[1.9], Δp=[0.2], Δs_a=[0.1; -0.1; 0.1], model_generator=create_jump_model_1), + "geq active bound change" => (p_a=[2.1], Δp=[-0.2], Δs_a=[-0.1; 0.1], model_generator=create_jump_model_2), + "leq active constraint change" => (p_a=[-1.9], Δp=[-0.2], Δs_a=[-0.1; 0.1; -0.1], model_generator=create_jump_model_3), + "leq active constraint change max" => (p_a=[1.9], Δp=[0.2], Δs_a=[0.1; -0.1; 0.1], model_generator=create_jump_model_4), + "geq active constraint change max" => (p_a=[1.9], Δp=[0.2], Δs_a=[0.1; -0.1; 0.1], model_generator=create_jump_model_5), +) + +function test_compute_derivatives_Analytical(DICT_PROBLEMS) + @testset "Compute Derivatives Analytical: $problem_name" for (problem_name, (p_a, Δp, Δs_a, model_generator)) in DICT_PROBLEMS + # OPT Problem + model, primal_vars, cons, params = model_generator() + eval_model_jump(model, primal_vars, cons, params, p_a) + # Compute derivatives + (Δs, sp_approx), evaluator, cons = compute_sensitivity(model, Δp; primal_vars, params) + # Check sensitivities + @test all(isapprox.(Δs[1:length(Δs_a)], Δs_a; atol = 1e-4)) + end +end + +################################################ +#= +# Test Sensitivity through finite differences +=# +################################################ + +function eval_model_jump(model, primal_vars, cons, params, p_val) + set_parameter_value.(params, p_val) + optimize!(model) + @assert is_solved_and_feasible(model) + return value.(primal_vars), dual.(cons), [dual.(LowerBoundRef(v)) for v in primal_vars if has_lower_bound(v)], [dual.(UpperBoundRef(v)) for v in primal_vars if has_upper_bound(v)] +end + +function stack_solution(cons, leq_locations, geq_locations, x, _λ, ν_L, ν_U) + ineq_locations = vcat(geq_locations, leq_locations) + return Float64[x; value.(get_slack_inequality.(cons[ineq_locations])); _λ; ν_L; _λ[geq_locations]; ν_U; _λ[leq_locations]] +end + +function print_wrong_sensitive(Δs, Δs_fd, primal_vars, cons, leq_locations, geq_locations) + ineq_locations = vcat(geq_locations, leq_locations) + println("Some sensitivities are not correct: \n") + # primal vars + num_primal_vars = length(primal_vars) + for (i, v) in enumerate(primal_vars) + if !isapprox(Δs[i], Δs_fd[i]; atol = 1e-6) + println("Primal var: ", v, " | Δs: ", Δs[i], " | Δs_fd: ", Δs_fd[i]) + end + end + # slack vars + num_slack_vars = length(ineq_locations) + num_w = num_slack_vars + num_primal_vars + for (i, c) in enumerate(cons[ineq_locations]) + if !isapprox(Δs[i + num_primal_vars], Δs_fd[i + num_primal_vars] ; atol = 1e-6) + println("Slack var: ", c, " | Δs: ", Δs[i + num_primal_vars], " | Δs_fd: ", Δs_fd[i + num_primal_vars]) + end + end + # dual vars + num_cons = length(cons) + for (i, c) in enumerate(cons) + if !isapprox(Δs[i + num_w], Δs_fd[i + num_w] ; atol = 1e-6) + println("Dual var: ", c, " | Δs: ", Δs[i + num_w], " | Δs_fd: ", Δs_fd[i + num_w]) + end + end + # dual lower bound primal vars + var_lower = [v for v in primal_vars if has_lower_bound(v)] + num_lower_bounds = length(var_lower) + for (i, v) in enumerate(var_lower) + if !isapprox(Δs[i + num_w + num_cons], Δs_fd[i + num_w + num_cons] ; atol = 1e-6) + lower_bound_ref = LowerBoundRef(v) + println("lower bound dual: ", lower_bound_ref, " | Δs: ", Δs[i + num_w + num_cons], " | Δs_fd: ", Δs_fd[i + num_w + num_cons]) + end + end + # dual lower bound slack vars + for (i, c) in enumerate(cons[geq_locations]) + if !isapprox(Δs[i + num_w + num_cons + num_lower_bounds], Δs_fd[i + num_w + num_cons + num_lower_bounds] ; atol = 1e-6) + println("lower bound slack dual: ", c, " | Δs: ", Δs[i + num_w + num_cons + num_lower_bounds], " | Δs_fd: ", Δs_fd[i + num_w + num_cons + num_lower_bounds]) + end + end + for (i, c) in enumerate(cons[leq_locations]) + if !isapprox(Δs[i + num_w + num_cons + num_lower_bounds + length(geq_locations)], Δs_fd[i + num_w + num_cons + num_lower_bounds + length(geq_locations)] ; atol = 1e-6) + println("upper bound slack dual: ", c, " | Δs: ", Δs[i + num_w + num_cons + num_lower_bounds + length(geq_locations)], " | Δs_fd: ", Δs_fd[i + num_w + num_cons + num_lower_bounds + length(geq_locations)]) + end + end + # dual upper bound primal vars + var_upper = [v for v in primal_vars if has_upper_bound(v)] + for (i, v) in enumerate(var_upper) + if !isapprox(Δs[i + num_w + num_cons + num_lower_bounds + num_slack_vars], Δs_fd[i + num_w + num_cons + num_lower_bounds + num_slack_vars] ; atol = 1e-6) + upper_bound_ref = UpperBoundRef(v) + println("upper bound dual: ", upper_bound_ref, " | Δs: ", Δs[i + num_w + num_cons + num_lower_bounds + num_slack_vars], " | Δs_fd: ", Δs_fd[i + num_w + num_cons + num_lower_bounds + num_slack_vars]) + end + end +end + +DICT_PROBLEMS_no_cc = Dict( + "QP_sIpopt" => (p_a=[4.5; 1.0], Δp=[0.001; 0.0], model_generator=create_nonlinear_jump_model_sipopt), + "NLP_1" => (p_a=[3.0; 2.0; 200], Δp=[0.001; 0.0; 0.0], model_generator=create_nonlinear_jump_model_1), + "NLP_1_2" => (p_a=[3.0; 2.0; 200], Δp=[0.0; 0.001; 0.0], model_generator=create_nonlinear_jump_model_1), + "NLP_1_3" => (p_a=[3.0; 2.0; 200], Δp=[0.0; 0.0; 0.001], model_generator=create_nonlinear_jump_model_1), + "NLP_1_4" => (p_a=[3.0; 2.0; 200], Δp=[0.1; 0.5; 0.5], model_generator=create_nonlinear_jump_model_1), + "NLP_1_4" => (p_a=[3.0; 2.0; 200], Δp=[0.5; -0.5; 0.1], model_generator=create_nonlinear_jump_model_1), + "NLP_2" => (p_a=[3.0; 2.0; 10], Δp=[0.01; 0.0; 0.0], model_generator=create_nonlinear_jump_model_2), + "NLP_2_2" => (p_a=[3.0; 2.0; 10], Δp=[-0.1; 0.0; 0.0], model_generator=create_nonlinear_jump_model_2), + "NLP_3" => (p_a=[3.0; 2.0; 10], Δp=[0.001; 0.0; 0.0], model_generator=create_nonlinear_jump_model_3), + "NLP_3_2" => (p_a=[3.0; 2.0; 10], Δp=[0.0; 0.001; 0.0], model_generator=create_nonlinear_jump_model_3), + "NLP_3_3" => (p_a=[3.0; 2.0; 10], Δp=[0.0; 0.0; 0.001], model_generator=create_nonlinear_jump_model_3), + "NLP_3_4" => (p_a=[3.0; 2.0; 10], Δp=[0.5; 0.001; 0.5], model_generator=create_nonlinear_jump_model_3), + "NLP_3_5" => (p_a=[3.0; 2.0; 10], Δp=[0.1; 0.3; 0.1], model_generator=create_nonlinear_jump_model_3), + "NLP_3_6" => (p_a=[3.0; 2.0; 10], Δp=[0.1; 0.2; -0.5], model_generator=create_nonlinear_jump_model_3), + "NLP_4" => (p_a=[1.0; 2.0; 100], Δp=[0.001; 0.0; 0.0], model_generator=create_nonlinear_jump_model_4), + "NLP_5" => (p_a=[1.0; 2.0; 100], Δp=[0.0; 0.001; 0.0], model_generator=create_nonlinear_jump_model_5), + "NLP_6" => (p_a=[100.0; 200.0], Δp=[0.2; 0.5], model_generator=create_nonlinear_jump_model_6), +) + + +DICT_PROBLEMS_cc = Dict( + "QP_JuMP" => (p_a=[1.0; 2.0; 100.0], Δp=[-0.5; 0.5; 0.1], model_generator=create_nonlinear_jump_model), + "QP_sIpopt2" => (p_a=[5.0; 1.0], Δp=[-0.5; 0.0], model_generator=create_nonlinear_jump_model_sipopt), +) + +function test_compute_derivatives_Finite_Diff(DICT_PROBLEMS, iscc=false) + # @testset "Compute Derivatives: $problem_name" + for (problem_name, (p_a, Δp, model_generator)) in DICT_PROBLEMS, ismin in [true, false] + # OPT Problem + model, primal_vars, cons, params = model_generator(;ismin=ismin) + eval_model_jump(model, primal_vars, cons, params, p_a) + println("$problem_name: ", model) + # Compute derivatives + # Δp = [0.001; 0.0; 0.0] + p_b = p_a .+ Δp + (Δs, sp_approx), evaluator, cons = compute_sensitivity(model, Δp; primal_vars, params) + leq_locations, geq_locations = find_inequealities(cons) + sa = stack_solution(cons, leq_locations, geq_locations, eval_model_jump(model, primal_vars, cons, params, p_a)...) + # Check derivatives using finite differences + ∂s_fd = FiniteDiff.finite_difference_jacobian((p) -> stack_solution(cons, leq_locations, geq_locations, eval_model_jump(model, primal_vars, cons, params, p)...), p_a) + Δs_fd = ∂s_fd * Δp + # actual solution + sp = stack_solution(cons, leq_locations, geq_locations, eval_model_jump(model, primal_vars, cons, params, p_b)...) + # Check sensitivities + num_important = length(primal_vars) + length(cons) + test_derivatives = all(isapprox.(Δs, Δs_fd; rtol = 1e-5, atol=1e-6)) + test_approx = all(isapprox.(sp[1:num_important], sp_approx[1:num_important]; rtol = 1e-5, atol=1e-6)) + if test_derivatives || (iscc && test_approx) + println("All sensitivities are correct") + elseif iscc && !test_approx + @show Δp + println("Fail Approximations") + print_wrong_sensitive(Δs, sp.-sa, primal_vars, cons, leq_locations, geq_locations) + else + @show Δp + print_wrong_sensitive(Δs, Δs_fd, primal_vars, cons, leq_locations, geq_locations) + end + println("--------------------") + end +end \ No newline at end of file diff --git a/opf.jl b/opf.jl new file mode 100644 index 00000000..e0d2697e --- /dev/null +++ b/opf.jl @@ -0,0 +1,433 @@ +#!/usr/bin/env julia +###### AC-OPF using JuMP ###### +# +# implementation reference: https://github.com/lanl-ansi/PowerModelsAnnex.jl/blob/master/src/model/ac-opf.jl +# only the built-in AD library is supported +# + +using PowerModels +using PGLib +import Ipopt +import JuMP +using LinearAlgebra +using ForwardDiff + +function build_gen(gen_bus, id, pmax, qmax) + return Dict{String, Any}( + "pg" => 0.5, + "model" => 2, + "shutdown" => 0.0, + "startup" => 0.0, + "qg" => 0.0, + "gen_bus" => gen_bus, + "pmax" => pmax, + "vg" => 1.0, + "mbase" => 100.0, + "source_id" => Any["gen", id], + "index" => id, + "cost" => [0.0, 0.0, 0.0], + "qmax" => qmax, + "gen_status" => 1, + "qmin" => - qmax, + "pmin" => 0.0, + "ncost" => 3, + ) +end + +function build_opf_model(data; add_param_load=false, solver=Ipopt.Optimizer) + # create ref + PowerModels.standardize_cost_terms!(data, order=2) + PowerModels.calc_thermal_limits!(data) + ref = PowerModels.build_ref(data)[:it][:pm][:nw][0] + + # create model + model = JuMP.Model(solver) + + JuMP.@variable(model, va[i in keys(ref[:bus])]) + JuMP.@variable(model, ref[:bus][i]["vmin"] <= vm[i in keys(ref[:bus])] <= ref[:bus][i]["vmax"], start=1.0) + for i in keys(ref[:bus]) + JuMP.set_name(va[i], "va_$i") + JuMP.set_name(vm[i], "vm_$i") + end + JuMP.@variable(model, ref[:gen][i]["pmin"] <= pg[i in keys(ref[:gen])] <= ref[:gen][i]["pmax"]) + JuMP.@variable(model, ref[:gen][i]["qmin"] <= qg[i in keys(ref[:gen])] <= ref[:gen][i]["qmax"]) + + for i in keys(ref[:gen]) + JuMP.set_name(pg[i], "pg_$i") + JuMP.set_name(qg[i], "qg_$i") + end + JuMP.@variable(model, -ref[:branch][l]["rate_a"] <= p[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"]) + JuMP.@variable(model, -ref[:branch][l]["rate_a"] <= q[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"]) + for (l, i, j) in ref[:arcs] + JuMP.set_name(p[(l, i, j)], "p_$(l)_$(i)_$(j)") + JuMP.set_name(q[(l, i, j)], "q_$(l)_$(i)_$(j)") + end + JuMP.@objective(model, Min, sum(gen["cost"][1]*pg[i]^2 + gen["cost"][2]*pg[i] + gen["cost"][3] for (i,gen) in ref[:gen])) + + for (i,bus) in ref[:ref_buses] + JuMP.@constraint(model, va[i] == 0) + end + + if add_param_load + @variable(model, pload[i in keys(ref[:bus])] ∈ MOI.Parameter.(0.0)) + @variable(model, qload[i in keys(ref[:bus])] ∈ MOI.Parameter.(0.0)) + else + pload = zeros(length(ref[:bus])) + qload = zeros(length(ref[:bus])) + end + + demand_equilibrium = Dict{Int, JuMP.ConstraintRef}() + for (i,bus) in ref[:bus] + bus_loads = [ref[:load][l] for l in ref[:bus_loads][i]] + bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][i]] + + active_balance = JuMP.@constraint(model, + sum(p[a] for a in ref[:bus_arcs][i]) == + sum(pg[g] for g in ref[:bus_gens][i]) - + sum(load["pd"] for load in bus_loads) - + pload[i] - + sum(shunt["gs"] for shunt in bus_shunts)*vm[i]^2 + ) + + JuMP.@constraint(model, + sum(q[a] for a in ref[:bus_arcs][i]) == + sum(qg[g] for g in ref[:bus_gens][i]) - + sum(load["qd"] for load in bus_loads) - + qload[i] + + sum(shunt["bs"] for shunt in bus_shunts)*vm[i]^2 + ) + + demand_equilibrium[i] = active_balance + end + + # Branch power flow physics and limit constraints + for (i,branch) in ref[:branch] + f_idx = (i, branch["f_bus"], branch["t_bus"]) + t_idx = (i, branch["t_bus"], branch["f_bus"]) + + p_fr = p[f_idx] + q_fr = q[f_idx] + p_to = p[t_idx] + q_to = q[t_idx] + + vm_fr = vm[branch["f_bus"]] + vm_to = vm[branch["t_bus"]] + va_fr = va[branch["f_bus"]] + va_to = va[branch["t_bus"]] + + g, b = PowerModels.calc_branch_y(branch) + tr, ti = PowerModels.calc_branch_t(branch) + ttm = tr^2 + ti^2 + g_fr = branch["g_fr"] + b_fr = branch["b_fr"] + g_to = branch["g_to"] + b_to = branch["b_to"] + + # From side of the branch flow + JuMP.@constraint(model, p_fr == (g+g_fr)/ttm*vm_fr^2 + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) ) + JuMP.@constraint(model, q_fr == -(b+b_fr)/ttm*vm_fr^2 - (-b*tr-g*ti)/ttm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/ttm*(vm_fr*vm_to*sin(va_fr-va_to)) ) + + # To side of the branch flow + JuMP.@constraint(model, p_to == (g+g_to)*vm_to^2 + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) ) + JuMP.@constraint(model, q_to == -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/ttm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-g*tr-b*ti)/ttm*(vm_to*vm_fr*sin(va_to-va_fr)) ) + + # Voltage angle difference limit + JuMP.@constraint(model, branch["angmin"] <= va_fr - va_to <= branch["angmax"]) + + # Apparent power limit, from side and to side + JuMP.@constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2) + JuMP.@constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2) + end + + return model, ref, demand_equilibrium, p, q, pg, qg, va, vm, pload, qload +end + +function build_bidding_opf_model(data; percen_bidding_nodes=0.1, solver=Ipopt.Optimizer) + data["basic_network"] = false + # add bidding generators + num_bidding_nodes = ceil(Int, length(data["bus"]) * percen_bidding_nodes) + bidding_nodes = parse.(Int, rand(keys(data["bus"]), num_bidding_nodes)) + existing_gens = maximum(parse.(Int, collect(keys(data["gen"])))) + pmax = maximum([data["gen"][g]["pmax"] for g in keys(data["gen"])]) + qmax = maximum([data["gen"][g]["qmax"] for g in keys(data["gen"])]) + bidding_gen_ids = existing_gens + 1:existing_gens + num_bidding_nodes + for (i, node) in enumerate(bidding_nodes) + data["gen"]["$(bidding_gen_ids[i])"] = build_gen(node, bidding_gen_ids[i], pmax, qmax) + end + data = make_basic_network(data) + total_market = sum([data["gen"][g]["pmax"] for g in keys(data["gen"])]) + + model, ref, demand_equilibrium, p, q, pg, qg, va, vm, pload, qload = build_opf_model(data; solver=solver) + + # add bids + JuMP.@variable(model, pg_bid[i in bidding_gen_ids] ∈ MOI.Parameter.(pmax * 0.5)) + @constraint(model, bid[i in bidding_gen_ids], pg[i] <= pg_bid[i]) + + model_variables = JuMP.num_variables(model) + + # for consistency with other solvers, skip the variable bounds in the constraint count + model_constraints = JuMP.num_constraints(model; count_variable_in_set_constraints = false) + + println("") + println("\033[1mSummary\033[0m") + # println(" case........: $(case_name)") + println(" variables...: $(model_variables)") + println(" constraints.: $(model_constraints)") + + println("") + + bidding_generators_dispatch = [pg[i] for i in bidding_gen_ids] + _pg_bid = [pg_bid[i] for i in bidding_gen_ids] + all_primal_variables = all_variables(model) + remaining_vars = setdiff(all_primal_variables, bidding_generators_dispatch) + remaining_vars = setdiff(remaining_vars, _pg_bid) + return Dict( + # "case" => case_name, + "model_variables" => remaining_vars, + "bidding_generators_dispatch" => bidding_generators_dispatch, + "bidding_lmps" => [demand_equilibrium[i] for i in bidding_nodes], + "model" => model, + "bid" => _pg_bid, + "pmax" => pmax, + "total_market" => total_market + ) +end + +function build_bidding_upper(num_bidding_nodes, pmax; solver=Ipopt.Optimizer) + upper_model = Model(solver) + + @variable(upper_model, 0 <= pg_bid[i=1:num_bidding_nodes] <= pmax) + @variable(upper_model, pg[i=1:num_bidding_nodes]) + @variable(upper_model, lambda[i=1:num_bidding_nodes]) + @objective(upper_model, Max, dot(lambda, pg)) + return upper_model, lambda, pg_bid, pg +end + +function memoize(foo::Function, n_outputs::Int) + last_x, last_f = nothing, nothing + function foo_i(i, x...) + if x !== last_x + ret = foo(x...) + if !isnothing(ret) + last_x, last_f = x, ret + end + end + return last_f[i] + end + return [(x...) -> foo_i(i, x...) for i in 1:n_outputs] +end + +function memoize!(foo::Function, n_outputs::Int) + last_x, last_f = nothing, nothing + function foo_i(i, g, x...) + if x !== last_x + ret = foo(g, x...) + if !isnothing(ret) + last_x, last_f = x, ret + end + end + g .= last_f[i] + return last_f[i] + end + return [(g, x...) -> foo_i(i, g, x...) for i in 1:n_outputs] +end + +function fdiff_derivatives(f::Function) + function ∇f(g::AbstractVector{T}, x::Vararg{T,N}) where {T,N} + FiniteDiff.finite_difference_gradient!(g, y -> f(y...), collect(x)) + return + end + function ∇²f(H::AbstractMatrix{T}, x::Vararg{T,N}) where {T,N} + h = FiniteDiff.finite_difference_hessian!(y -> f(y...), collect(x)) + for i in 1:N, j in 1:i + H[i, j] = h[i, j] + end + return + end + return ∇f, ∇²f +end + +function test_bilevel_ac_strategic_bidding(data; percen_bidding_nodes=0.1, Δp=0.0001, solver_upper=Ipopt.Optimizer, solver_lower=Ipopt.Optimizer) + # test derivative of the dual of the demand equilibrium constraint + data = build_bidding_opf_model(data; percen_bidding_nodes=percen_bidding_nodes, solver=solver_lower) + pmax = data["pmax"] + primal_vars = [data["bidding_generators_dispatch"]; data["model_variables"]] + num_bidding_nodes = length(data["bidding_generators_dispatch"]) + set_parameter_value.(data["bid"], 0.01) + JuMP.optimize!(data["model"]) + evaluator, cons = create_evaluator(data["model"]; x=[primal_vars; data["bid"]]) + leq_locations, geq_locations = find_inequealities(cons) + num_ineq = length(leq_locations) + length(geq_locations) + num_primal = length(primal_vars) + bidding_lmps_index = [findall(x -> x == i, cons)[1] for i in data["bidding_lmps"]] + if !isnothing(Δp) + Δp = fill(Δp, num_bidding_nodes) + end + # Δp = rand(-pmax*0.1:0.001:pmax*0.1, num_bidding_nodes) + # Δs, sp = compute_sensitivity(evaluator, cons, Δp; primal_vars=primal_vars, params=data["bid"]) + + # set_parameter_value.(data["bid"], 0.5 .+ Δp) + # JuMP.optimize!(data["model"]) + + # @test dual.(data["bidding_lmps"]) ≈ Δs[(num_primal + num_ineq) .+ bidding_lmps_index] + + # test bilevel strategic bidding + upper_model, lambda, pg_bid, pg = build_bidding_upper(num_bidding_nodes, pmax; solver=solver_upper) + evaluator, cons = create_evaluator(data["model"]; x=[primal_vars; data["bid"]]) + + function f(pg_bid_val...) + set_parameter_value.(data["bid"], pg_bid_val) + JuMP.optimize!(data["model"]) + if !is_solved_and_feasible(data["model"]) + return nothing + end + return [value.(data["bidding_generators_dispatch"]); -dual.(data["bidding_lmps"])] + end + + function ∇f(g::AbstractVector{T}, pg_bid_val...) where {T} + if !all(value.(data["bid"]) .== pg_bid_val) + set_parameter_value.(data["bid"], pg_bid_val) + JuMP.optimize!(data["model"]) + @assert is_solved_and_feasible(data["model"]) + end + + Δs, sp = compute_sensitivity(evaluator, cons, Δp; primal_vars=primal_vars, params=data["bid"]) + if isnothing(Δp) + Δs = Δs * ones(size(Δs, 2)) + end + for i in 1:num_bidding_nodes + g[i] = Δs[i] + end + for (i, b_idx) in enumerate(bidding_lmps_index) + g[i] = -Δs[num_primal + num_ineq + b_idx] + end + return g + # return [Δs[1:num_bidding_nodes]; -(Δs[(num_primal + num_ineq) .+ bidding_lmps_index])] + end + + memoized_f = memoize(f, 2 * num_bidding_nodes) + memoized_∇f = !isnothing(Δp) && iszero(Δp) ? [(args...) -> nothing for i in 1:2 * num_bidding_nodes] : memoize!(∇f, 2 * num_bidding_nodes) + + for i in 1:num_bidding_nodes + op_pg = add_nonlinear_operator(upper_model, num_bidding_nodes, memoized_f[i], memoized_∇f[i]; name = Symbol("op_pg_$i")) + # op_pg = add_nonlinear_operator(upper_model, num_bidding_nodes, memoized_f[i], fdiff_derivatives(memoized_f[i])...; name = Symbol("op_pg_$i")) + @constraint(upper_model, pg[i] == op_pg(pg_bid...)) + end + + for i in 1:num_bidding_nodes + op_lambda = add_nonlinear_operator(upper_model, num_bidding_nodes, memoized_f[num_bidding_nodes + i], memoized_∇f[num_bidding_nodes + i]; name = Symbol("op_lambda_$i")) + # op_lambda = add_nonlinear_operator(upper_model, num_bidding_nodes, memoized_f[num_bidding_nodes + i], fdiff_derivatives(memoized_f[num_bidding_nodes + i])...; name = Symbol("op_lambda_$i")) + @constraint(upper_model, lambda[i] == op_lambda(pg_bid...)) + end + + JuMP.optimize!(upper_model) + + println("Status: ", termination_status(upper_model)) + println("Objective: ", objective_value(upper_model)) + println("Duals: ", value.(lambda)) + println("Bids: ", value.(pg_bid)) + println("Dispatch: ", value.(pg)) +end + +function test_bidding_nlopt(data; percen_bidding_nodes=0.1, Δp=0.0001, solver_upper=:LD_MMA, solver_lower=optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0), + max_eval=100, pmax_multiplier=1.0, bool_trace=false +) + data = build_bidding_opf_model(data; percen_bidding_nodes=percen_bidding_nodes, solver=solver_lower) + pmax = data["pmax"] * pmax_multiplier + primal_vars = [data["bidding_generators_dispatch"]; data["model_variables"]] + num_bidding_nodes = length(data["bidding_generators_dispatch"]) + set_parameter_value.(data["bid"], 0.01) + JuMP.optimize!(data["model"]) + evaluator, cons = create_evaluator(data["model"]; x=[primal_vars; data["bid"]]) + leq_locations, geq_locations = find_inequealities(cons) + num_ineq = length(leq_locations) + length(geq_locations) + num_primal = length(primal_vars) + bidding_lmps_index = [findall(x -> x == i, cons)[1] for i in data["bidding_lmps"]] + if !isnothing(Δp) + Δp = fill(Δp, num_bidding_nodes) + end + + evaluator, cons = create_evaluator(data["model"]; x=[primal_vars; data["bid"]]) + + function f(pg_bid_val...) + set_parameter_value.(data["bid"], pg_bid_val) + JuMP.optimize!(data["model"]) + @assert is_solved_and_feasible(data["model"]) + + return value.(data["bidding_generators_dispatch"]), -dual.(data["bidding_lmps"]) + end + + function ∇f(pg_bid_val...) + @assert all(value.(data["bid"]) .== pg_bid_val) + + Δs, sp = compute_sensitivity(evaluator, cons, Δp; primal_vars=primal_vars, params=data["bid"]) + if isnothing(Δp) + Δs = Δs * ones(size(Δs, 2)) + end + + return Δs[1:num_bidding_nodes], -(Δs[(num_primal + num_ineq) .+ bidding_lmps_index]) + end + + if !isnothing(Δp) && iszero(Δp) + _∇f = (pg_bid_val...) -> (zeros(num_bidding_nodes), zeros(length(bidding_lmps_index))) + else + _∇f = ∇f + end + + trace = Any[] + function my_objective_fn(pg_bid_val::Vector, grad::Vector) + pg, lmps = f(pg_bid_val...) + Δpg, Δlmps = _∇f(pg_bid_val...) + value = dot(lmps, pg) + if length(grad) > 0 + grad .= dot(Δlmps, pg) .+ dot(lmps, Δpg) + end + # println("Objective: ", value) + if bool_trace + push!(trace, sum(pg_bid_val) => value) + end + return value + end + + opt = NLopt.Opt(solver_upper, num_bidding_nodes) + if solver_upper == :G_MLSL_LDS + NLopt.nlopt_set_local_optimizer(opt, NLopt.Opt(:LN_NELDERMEAD, num_bidding_nodes)) + end + NLopt.lower_bounds!(opt, fill(0.0, num_bidding_nodes)) + NLopt.upper_bounds!(opt, fill(pmax, num_bidding_nodes)) + NLopt.xtol_rel!(opt, 1e-4) + maxeval!(opt, max_eval) + NLopt.max_objective!(opt, my_objective_fn) + max_f, opt_x, ret = NLopt.optimize(opt, rand(num_bidding_nodes)) + num_evals = NLopt.numevals(opt) + println("Status: ", ret) + println("Objective: ", max_f) + println("Bids: ", opt_x[1:num_bidding_nodes]) + println("Duals: ", opt_x[num_bidding_nodes+1:end]) + println("Number of evaluations: ", num_evals) + return max_f, num_evals, trace, (sum(opt_x[1:num_bidding_nodes]) / data["total_market"]) * 100, ret +end + +function sesitivity_load(data; Δp=nothing, solver=Ipopt.Optimizer) + # data = make_basic_network(pglib(case_name)) + + model, ref, demand_equilibrium, p, q, pg, qg, va, vm, pload, qload = build_opf_model(data; solver=solver, add_param_load=true) + num_bus = length(ref[:bus]) + demand_equilibrium = [demand_equilibrium[i] for i in 1:num_bus] + + JuMP.optimize!(model) + + params = vcat(pload.data, qload.data) + all_primal_variables = all_variables(model) + primal_vars = setdiff(all_primal_variables, params) + + evaluator, cons = create_evaluator(model; x=[primal_vars; params]) + + leq_locations, geq_locations = find_inequealities(cons) + num_ineq = length(leq_locations) + length(geq_locations) + num_primal = length(primal_vars) + lmps_index = [findall(x -> x == i, cons)[1] for i in demand_equilibrium] + + Δs, sp = compute_sensitivity(evaluator, cons, Δp; primal_vars=primal_vars, params=params) + return Δs[1:num_primal, :], Δs[(num_primal + num_ineq) .+ lmps_index, :] +end diff --git a/problem_list.jl b/problem_list.jl new file mode 100644 index 00000000..ef551a32 --- /dev/null +++ b/problem_list.jl @@ -0,0 +1,340 @@ +using JuMP +using Ipopt + +################################################ +#= +From JuMP Tutorial for Querying Hessians: +https://github.com/jump-dev/JuMP.jl/blob/301d46e81cb66c74c6e22cd89fb89ced740f157b/docs/src/tutorials/nonlinear/querying_hessians.jl#L67-L72 +=# +################################################ +function create_nonlinear_jump_model(;ismin=true) + model = Model(Ipopt.Optimizer) + set_silent(model) + @variable(model, p ∈ MOI.Parameter(1.0)) + @variable(model, p2 ∈ MOI.Parameter(2.0)) + @variable(model, p3 ∈ MOI.Parameter(100.0)) + @variable(model, x[i = 1:2], start = -i) + @constraint(model, g_1, x[1]^2 <= p) + @constraint(model, g_2, p * (x[1] + x[2])^2 <= p2) + if ismin + @objective(model, Min, (1 - x[1])^2 + p3 * (x[2] - x[1]^2)^2) + else + @objective(model, Max, -(1 - x[1])^2 - p3 * (x[2] - x[1]^2)^2) + end + + return model, x, [g_1; g_2], [p; p2; p3] +end + + +################################################ +#= +From sIpopt paper: https://optimization-online.org/2011/04/3008/ +=# +################################################ + +function create_nonlinear_jump_model_sipopt(;ismin = true) + model = Model(Ipopt.Optimizer) + set_silent(model) + @variable(model, p1 ∈ MOI.Parameter(4.5)) + @variable(model, p2 ∈ MOI.Parameter(1.0)) + @variable(model, x[i = 1:3] >= 0, start = -i) + @constraint(model, g_1, 6 * x[1] + 3 * x[2] + 2 * x[3] - p1 == 0) + @constraint(model, g_2, p2 * x[1] + x[2] - x[3] - 1 == 0) + if ismin + @objective(model, Min, x[1]^2 + x[2]^2 + x[3]^2) + else + @objective(model, Max, -x[1]^2 - x[2]^2 - x[3]^2) + end + return model, x, [g_1; g_2], [p1; p2] +end + +################################################ +#= +Simple Problems +=# +################################################ + + +function create_jump_model_1(p_val = [1.5]) + model = Model(Ipopt.Optimizer) + set_silent(model) + + # Parameters + @variable(model, p ∈ MOI.Parameter(p_val[1])) + + # Variables + @variable(model, x) + + # Constraints + @constraint(model, con1, x >= p) + @constraint(model, con2, x >= 2) + @objective(model, Min, x^2) + + return model, [x], [con1; con2], [p] +end + +function create_jump_model_2(p_val = [1.5]) + model = Model(Ipopt.Optimizer) + set_silent(model) + + # Parameters + @variable(model, p ∈ MOI.Parameter(p_val[1])) + + # Variables + @variable(model, x >= 2.0) + + # Constraints + @constraint(model, con1, x >= p) + @objective(model, Min, x^2) + + return model, [x], [con1], [p] +end + +function create_jump_model_3(p_val = [-1.5]) + model = Model(Ipopt.Optimizer) + set_silent(model) + + # Parameters + @variable(model, p ∈ MOI.Parameter(p_val[1])) + + # Variables + @variable(model, x) + + # Constraints + @constraint(model, con1, x <= p) + @constraint(model, con2, x <= -2) + @objective(model, Min, -x) + + return model, [x], [con1; con2], [p] +end + +function create_jump_model_4(p_val = [1.5]) + model = Model(Ipopt.Optimizer) + set_silent(model) + + # Parameters + @variable(model, p ∈ MOI.Parameter(p_val[1])) + + # Variables + @variable(model, x) + + # Constraints + @constraint(model, con1, x <= p) + @constraint(model, con2, x <= 2) + @objective(model, Max, x) + + return model, [x], [con1; con2], [p] +end + +function create_jump_model_5(p_val = [1.5]) + model = Model(Ipopt.Optimizer) + set_silent(model) + + # Parameters + @variable(model, p ∈ MOI.Parameter(p_val[1])) + + # Variables + @variable(model, x) + + # Constraints + @constraint(model, con1, x >= p) + @constraint(model, con2, x >= 2) + @objective(model, Max, -x) + + return model, [x], [con1; con2], [p] +end + +# Softmax model +h(y) = - sum(y .* log.(y)) +softmax(x) = exp.(x) / sum(exp.(x)) +function create_jump_model_6(p_a = collect(1.0:0.1:2.0)) + model = Model(Ipopt.Optimizer) + set_silent(model) + + # Parameters + @variable(model, x[i=1:length(p_a)] ∈ MOI.Parameter.(p_a)) + + # Variables + @variable(model, y[1:length(p_a)] >= 0.0) + + # Constraints + @constraint(model, con1, sum(y) == 1) + @constraint(model, con2[i=1:length(x)], y[i] <= 1) + + # Objective + @objective(model, Max, dot(x, y) + h(y)) + + return model, y, [con1; con2], x +end + +function create_jump_model_7(p_val = [1.5], g = sin) + model = Model(Ipopt.Optimizer) + set_silent(model) + + # Parameters + @variable(model, p ∈ MOI.Parameter(p_val[1])) + + # Variables + @variable(model, x) + + # Constraints + @constraint(model, con1, x * g(p) == 1) + @objective(model, Min, 0) + + return model, [x], [con1], [p] +end + +################################################ +#= +Non Linear Problems +=# +################################################ + + +function create_nonlinear_jump_model_1(p_val = [1.0; 2.0; 100]; ismin = true) + model = Model(Ipopt.Optimizer) + set_silent(model) + + # Parameters + @variable(model, p[i=1:3] ∈ MOI.Parameter.(p_val)) + + # Variables + @variable(model, x) + @variable(model, y) + + # Constraints + @constraint(model, con1, y >= p[1]*sin(x)) # NLP Constraint + @constraint(model, con2, x + y == p[1]) + @constraint(model, con3, p[2] * x >= 0.1) + if ismin + @objective(model, Min, (1 - x)^2 + p[3] * (y - x^2)^2) # NLP Objective + else + @objective(model, Max, -(1 - x)^2 - p[3] * (y - x^2)^2) # NLP Objective + end + + return model, [x; y], [con1; con2; con3], p +end + +function create_nonlinear_jump_model_2(p_val = [3.0; 2.0; 10]; ismin = true) + model = Model(Ipopt.Optimizer) + set_silent(model) + + # Parameters + @variable(model, p[i=1:3] ∈ MOI.Parameter.(p_val)) + + # Variables + @variable(model, x <= 10) + @variable(model, y) + + # Constraints + @constraint(model, con1, y >= p[1]*sin(x)) # NLP Constraint + @constraint(model, con2, x + y == p[1]) + @constraint(model, con3, p[2] * x >= 0.1) + if ismin + @objective(model, Min, (1 - x)^2 + p[3] * (y - x^2)^2) # NLP Objective + else + @objective(model, Max, -(1 - x)^2 - p[3] * (y - x^2)^2) # NLP Objective + end + + return model, [x; y], [con1; con2; con3], p +end + +function create_nonlinear_jump_model_3(p_val = [3.0; 2.0; 10]; ismin = true) + model = Model(Ipopt.Optimizer) + set_silent(model) + + # Parameters + @variable(model, p[i=1:3] ∈ MOI.Parameter.(p_val)) + + # Variables + @variable(model, x <= 10) + @variable(model, y) + + # Constraints + @constraint(model, con1, y >= p[1]*sin(x)) # NLP Constraint + @constraint(model, con2, x + y == p[1]) + @constraint(model, con3, p[2] * x >= 0.1) + if ismin + @objective(model, Min, (1 - x)^2 + p[3] * (y - x^2)^2) # NLP Objective + else + @objective(model, Max, -(1 - x)^2 - p[3] * (y - x^2)^2) # NLP Objective + end + return model, [x; y], [con1; con2; con3], p +end + +function create_nonlinear_jump_model_4(p_val = [1.0; 2.0; 100]; ismin = true) + model = Model(Ipopt.Optimizer) + set_silent(model) + + # Parameters + @variable(model, p[i=1:3] ∈ MOI.Parameter.(p_val)) + + # Variables + @variable(model, x) + @variable(model, y) + + # Constraints + @constraint(model, con0, x == p[1] - 0.5) + @constraint(model, con1, y >= p[1]*sin(x)) # NLP Constraint + @constraint(model, con2, x + y == p[1]) + @constraint(model, con3, p[2] * x >= 0.1) + if ismin + @objective(model, Min, (1 - x)^2 + p[3] * (y - x^2)^2) # NLP Objective + else + @objective(model, Max, -(1 - x)^2 - p[3] * (y - x^2)^2) # NLP Objective + end + + return model, [x; y], [con1; con2; con3], p +end + +function create_nonlinear_jump_model_5(p_val = [1.0; 2.0; 100]; ismin = true) + model = Model(Ipopt.Optimizer) + set_silent(model) + + # Parameters + @variable(model, p[i=1:3] ∈ MOI.Parameter.(p_val)) + + # Variables + @variable(model, x) + @variable(model, y) + + # Constraints + fix(x, 0.5) + con0 = JuMP.FixRef(x) + @constraint(model, con1, y >= p[1]*sin(x)) # NLP Constraint + @constraint(model, con2, x + y == p[1]) + @constraint(model, con3, p[2] * x >= 0.1) + if ismin + @objective(model, Min, (1 - x)^2 + p[3] * (y - x^2)^2) # NLP Objective + else + @objective(model, Max, -(1 - x)^2 - p[3] * (y - x^2)^2) # NLP Objective + end + + return model, [x; y], [con0; con1; con2; con3], p +end + +function create_nonlinear_jump_model_6(p_val = [100.0; 200.0]; ismin = true) + model = Model(Ipopt.Optimizer) + set_silent(model) + + # Parameters + @variable(model, p[i=1:2] ∈ MOI.Parameter.(p_val)) + + # Variables + @variable(model, x[i=1:2]) + @variable(model, z) # >= 2.0) + @variable(model, w) # <= 3.0) + # @variable(model, f[1:2]) + + # Constraints + @constraint(model, con1, x[2] - 0.0001 * x[1]^2 - 0.2 * z^2 - 0.3 * w^2 >= p[1] + 1) + @constraint(model, con2, x[1] + 0.001 * x[2]^2 + 0.5 * w^2 + 0.4 * z^2 <= 10 * p[1] + 2) + @constraint(model, con3, z^2 + w^2 == 13) + if ismin + @objective(model, Min, x[2] - x[1] + z - w) + else + @objective(model, Max, -x[2] + x[1] - z + w) + end + + return model, [x; z; w], [con2; con3], p +end diff --git a/results/complexity_mpec.csv b/results/complexity_mpec.csv new file mode 100644 index 00000000..07674d4b --- /dev/null +++ b/results/complexity_mpec.csv @@ -0,0 +1,12 @@ +nodes,nlp,exact +7,1.9666185609002593,2.7846410583761947 +12,2.0380107267563576,3.1072689575976105 +17,3.09061303995118,10.24712725907874 +22,2.297800196291057,13.512257289968211 +27,2.343125283792045,20.846918358124633 +32,356.02001668056715,154373.9398308114 +37,437.41794076766814,40769.94301293702 +42,276.7964300892478,452999.32044198894 +47,706.4255439924314,6.779457147587512e6 +52,745.8290398126463,3.3012040145018916e6 +57,670.5347723909827,3.824821345983114e6 diff --git a/results/complexity_mpec.pdf b/results/complexity_mpec.pdf new file mode 100644 index 00000000..885a022d Binary files /dev/null and b/results/complexity_mpec.pdf differ diff --git a/results/strategic_bidding_nlopt_pglib_opf_case1354_pegase.csv b/results/strategic_bidding_nlopt_pglib_opf_case1354_pegase.csv new file mode 100644 index 00000000..cf51f812 --- /dev/null +++ b/results/strategic_bidding_nlopt_pglib_opf_case1354_pegase.csv @@ -0,0 +1,81 @@ +solver_upper,solver_lower,Δp,seed,profit,market_share,num_evals,time,status +LD_MMA,Ipopt,nothing,1,379521.42897568096,4.463315414691349,54,1026.4853658676147,1 +LD_MMA,Ipopt,nothing,2,407060.93605103396,4.542909921306311,58,1076.0979299545288,1 +LD_MMA,Ipopt,nothing,3,410312.7304017537,5.129981875986086,61,1219.9480831623077,1 +LD_MMA,Ipopt,nothing,4,382648.40172933176,4.570043129002875,58,1093.451698064804,1 +LD_MMA,Ipopt,nothing,5,400129.3840004562,4.435130100704832,55,1004.400083065033,1 +LD_MMA,Ipopt,nothing,6,398863.17915605806,4.711472480818748,54,1025.8297719955444,1 +LD_MMA,Ipopt,nothing,7,366467.5619779321,4.418414297264759,61,1150.896108865738,1 +LD_MMA,Ipopt,nothing,8,404302.97117886203,4.759032820654829,50,931.877956867218,1 +LD_MMA,Ipopt,nothing,9,412951.31729381235,4.416180195739012,58,1071.8616619110107,1 +LD_MMA,Ipopt,nothing,10,391522.83322145656,4.46357445793139,60,1151.4680500030518,1 +LN_BOBYQA,Ipopt,0.0,1,168820.5939024172,1.004908393034264,100,523.141793012619,1 +LN_BOBYQA,Ipopt,0.0,2,158509.44145516364,0.9567626049305376,100,544.5915019512177,1 +LN_BOBYQA,Ipopt,0.0,3,166297.78580567916,1.0253714817670256,100,532.0673048496246,1 +LN_BOBYQA,Ipopt,0.0,4,158518.2600053562,0.9668237599289118,100,484.60249185562134,1 +LN_BOBYQA,Ipopt,0.0,5,162552.88193519256,0.981494279687558,100,535.4164979457855,1 +LN_BOBYQA,Ipopt,0.0,6,158328.14431752235,0.9888221309137759,100,497.9929940700531,1 +LN_BOBYQA,Ipopt,0.0,7,163347.59302177207,0.991401346808828,100,530.5398418903351,1 +LN_BOBYQA,Ipopt,0.0,8,162521.7493274639,0.9920052077878728,100,671.4308149814606,1 +LN_BOBYQA,Ipopt,0.0,9,163639.85254358838,0.9906836031732348,100,507.2492690086365,1 +LN_BOBYQA,Ipopt,0.0,10,157648.53408377702,0.9426881226241853,100,520.3809530735016,1 +LD_LBFGS,Ipopt,nothing,1,379533.928182014,4.463652068653421,71,1367.5183429718018,1 +LD_LBFGS,Ipopt,nothing,2,407216.4369955502,4.540208357985074,62,1169.3655450344086,1 +LD_LBFGS,Ipopt,nothing,3,410318.18203851796,5.130083101101336,63,1199.300220966339,1 +LD_LBFGS,Ipopt,nothing,4,375902.4215017463,4.661620288652698,27,538.6267030239105,1 +LD_LBFGS,Ipopt,nothing,5,404906.1092889891,4.088859405877921,19,351.5564959049225,1 +LD_LBFGS,Ipopt,nothing,6,402324.6784042329,4.668206955671465,49,909.8845860958099,1 +LD_LBFGS,Ipopt,nothing,7,406897.83261238714,4.206412789263682,49,929.1883580684662,1 +LD_LBFGS,Ipopt,nothing,8,404483.36856513243,4.75736385353077,51,971.0387299060822,1 +LD_LBFGS,Ipopt,nothing,9,431654.9285262448,4.334209147267336,67,1284.628026008606,1 +LD_LBFGS,Ipopt,nothing,10,401019.60731385625,4.229152525971723,28,532.8773579597473,1 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,1,379543.70165973814,4.463943219982909,107,2173.238525867462,4 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,6,402328.87271081784,4.668319190463052,85,1627.645066022873,4 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,8,406133.5530688092,4.3972176556105715,52,1025.2977900505066,4 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,9,431642.9372922031,4.333970990756671,50,1007.0831871032715,4 +LN_COBYLA,Ipopt,0.0,1,178359.74463298658,1.1200533758856948,100,593.1553659439087,5 +LN_COBYLA,Ipopt,0.0,2,212254.08379188183,1.4180001908515154,100,659.3406028747559,5 +LN_COBYLA,Ipopt,0.0,3,200066.04751296729,1.3332900353516668,100,595.0377469062805,5 +LN_COBYLA,Ipopt,0.0,4,218238.51484798797,1.4738619006311162,100,627.186784029007,5 +LN_COBYLA,Ipopt,0.0,5,203024.17273709437,1.3468078531035172,100,620.616494178772,5 +LN_COBYLA,Ipopt,0.0,6,191606.14917259244,1.2721033726111335,100,619.6530439853668,5 +LN_COBYLA,Ipopt,0.0,7,216634.18097193688,1.4533418695121343,100,620.3641991615295,5 +LN_COBYLA,Ipopt,0.0,8,220015.7631276374,1.4754972444825114,100,645.8744649887085,5 +LN_COBYLA,Ipopt,0.0,9,214629.68380120033,1.4855060084960066,100,580.2151620388031,5 +LN_COBYLA,Ipopt,0.0,10,189925.5983051968,1.2210286163994604,100,572.2593560218811,5 +LD_CCSAQ,Ipopt,nothing,1,286001.84954590956,6.691799689272085,44,869.4242498874664,4 +LD_CCSAQ,Ipopt,nothing,2,288247.7411749264,6.3875171811775004,35,649.8028280735016,4 +LD_CCSAQ,Ipopt,nothing,3,310628.608337936,6.372146554628574,58,1256.3902189731598,4 +LD_CCSAQ,Ipopt,nothing,5,298608.34558092774,6.8261547881138975,53,1064.9097120761871,4 +LD_CCSAQ,Ipopt,nothing,6,275603.01813745056,6.3531033574974884,55,1096.9445140361786,4 +LD_CCSAQ,Ipopt,nothing,7,289939.19647312385,6.440085995867354,37,724.0602381229401,4 +LD_CCSAQ,Ipopt,nothing,9,297084.5813973348,6.798650024586453,57,1227.8032429218292,4 +LD_SLSQP,Ipopt,nothing,1,376542.9974470474,3.2596009023573256,100,2095.886658191681,5 +LD_SLSQP,Ipopt,nothing,2,156706.26214407926,0.952168590574678,1,29.89198112487793,4 +LD_SLSQP,Ipopt,nothing,3,394748.4905896696,4.871532266285069,100,2073.313261985779,5 +LD_SLSQP,Ipopt,nothing,4,156773.4713477757,0.9610650569992588,1,29.20994520187378,4 +LD_SLSQP,Ipopt,nothing,5,161022.97906556918,0.9814364867894274,1,29.023494005203247,4 +LD_SLSQP,Ipopt,nothing,6,156720.88459197825,0.9791306609134534,1,28.679763793945312,4 +LD_SLSQP,Ipopt,nothing,7,161525.52344752845,0.983549520184631,1,29.40539813041687,4 +LD_SLSQP,Ipopt,nothing,8,394271.15137901035,4.110846793336604,100,2074.1711959838867,5 +LD_SLSQP,Ipopt,nothing,9,161725.68930725995,0.9896961020228838,1,30.21243405342102,4 +LD_SLSQP,Ipopt,nothing,10,385325.4647475914,3.744578744748129,100,2041.07581615448,5 +LN_NELDERMEAD,Ipopt,0.0,1,178359.74343438007,1.1200533758856948,100,511.4894938468933,5 +LN_NELDERMEAD,Ipopt,0.0,2,212254.08379190153,1.4180001908515154,100,550.9402630329132,5 +LN_NELDERMEAD,Ipopt,0.0,3,200066.04751296242,1.3332900353516668,100,512.9112930297852,5 +LN_NELDERMEAD,Ipopt,0.0,4,218238.51484798538,1.4738619006311162,100,501.04356598854065,5 +LN_NELDERMEAD,Ipopt,0.0,5,203024.17273709504,1.3468078531035172,100,521.932450056076,5 +LN_NELDERMEAD,Ipopt,0.0,6,191606.14917260138,1.2721033726111335,100,524.1365430355072,5 +LN_NELDERMEAD,Ipopt,0.0,7,216634.18097194025,1.4533418695121343,100,515.4771430492401,5 +LN_NELDERMEAD,Ipopt,0.0,9,214629.68380120143,1.4855060084960066,100,501.47680020332336,5 +LN_NELDERMEAD,Ipopt,0.0,10,189925.59830520058,1.2210286163994604,100,501.95699095726013,5 +LN_NEWUOA_BOUND,Ipopt,0.0,1,167196.39043391688,1.0022938482589328,100,512.7025361061096,5 +LN_NEWUOA_BOUND,Ipopt,0.0,2,156722.70401364815,0.952168590574678,100,535.0520720481873,5 +LN_NEWUOA_BOUND,Ipopt,0.0,3,164443.08291949987,1.0203179439844057,100,493.93674206733704,5 +LN_NEWUOA_BOUND,Ipopt,0.0,4,156794.56570489882,0.9610650569992588,100,482.8457419872284,5 +LN_NEWUOA_BOUND,Ipopt,0.0,5,161033.91621436726,0.9814364867894274,100,534.1026439666748,5 +LN_NEWUOA_BOUND,Ipopt,0.0,6,156725.4964958931,0.9791306609134534,100,475.3539800643921,5 +LN_NEWUOA_BOUND,Ipopt,0.0,7,161535.26400990074,0.983549520184631,100,499.9822299480438,5 +LN_NEWUOA_BOUND,Ipopt,0.0,8,160979.44597772305,0.9917204762644776,100,574.3065629005432,5 +LN_NEWUOA_BOUND,Ipopt,0.0,9,161764.20237385298,0.9896961020228838,100,493.1457760334015,5 +LN_NEWUOA_BOUND,Ipopt,0.0,10,155868.95779836978,0.942203145482624,100,510.560693025589,5 diff --git a/results/strategic_bidding_nlopt_pglib_opf_case1354_pegase.pdf b/results/strategic_bidding_nlopt_pglib_opf_case1354_pegase.pdf new file mode 100644 index 00000000..9f6e04ea Binary files /dev/null and b/results/strategic_bidding_nlopt_pglib_opf_case1354_pegase.pdf differ diff --git a/results/strategic_bidding_nlopt_pglib_opf_case2000_goc.csv b/results/strategic_bidding_nlopt_pglib_opf_case2000_goc.csv new file mode 100644 index 00000000..2dffa9fe --- /dev/null +++ b/results/strategic_bidding_nlopt_pglib_opf_case2000_goc.csv @@ -0,0 +1,73 @@ +solver_upper,solver_lower,Δp,seed,profit,market_share,num_evals,time,status +LD_MMA,Ipopt,nothing,1,369523.35100522,6.281331171600099,57,2927.388349056244,4 +LD_MMA,Ipopt,nothing,2,319324.6427443133,6.434705664578061,57,2780.2165780067444,4 +LD_MMA,Ipopt,nothing,3,360352.16414853133,6.5808096766567,60,3062.843197107315,4 +LD_MMA,Ipopt,nothing,4,363082.0923797429,6.460565006687831,60,3123.9112889766693,4 +LD_MMA,Ipopt,nothing,5,374396.0648059043,6.328558034627275,57,2961.078145980835,4 +LD_MMA,Ipopt,nothing,6,358114.98152800533,6.191722678203054,58,2993.763321876526,4 +LD_MMA,Ipopt,nothing,7,358894.43974197324,6.275367062460499,58,2930.235079050064,4 +LD_MMA,Ipopt,nothing,8,346924.9173690083,6.330849145471191,61,3128.324747800827,4 +LD_MMA,Ipopt,nothing,9,353932.7637797668,6.591812236705967,58,2978.2513270378113,4 +LD_MMA,Ipopt,nothing,10,360641.74436154397,6.394534482828754,58,2845.8516099452972,4 +LN_BOBYQA,Ipopt,0.0,1,297977.62147699925,3.8930735360049527,100,1144.0847730636597,5 +LN_BOBYQA,Ipopt,0.0,2,299355.77298776305,3.9404747232021657,100,1082.5064489841461,5 +LN_BOBYQA,Ipopt,0.0,3,277492.64827409753,3.759131393699533,100,1122.313206911087,5 +LN_BOBYQA,Ipopt,0.0,4,307781.33752304706,4.059088566045905,100,1152.3991768360138,5 +LN_BOBYQA,Ipopt,0.0,5,295187.4306268543,3.859893460772294,100,1182.9611599445343,5 +LN_BOBYQA,Ipopt,0.0,6,273338.4756316369,3.6000920228228708,100,954.2082068920135,5 +LN_BOBYQA,Ipopt,0.0,7,295616.27458341373,3.7874468131431374,100,999.1024839878082,5 +LN_BOBYQA,Ipopt,0.0,8,301629.43296452175,3.904028118971772,100,973.1779978275299,5 +LN_BOBYQA,Ipopt,0.0,9,286589.33802552376,3.8070028445188564,100,993.2643630504608,5 +LN_BOBYQA,Ipopt,0.0,10,306357.57064501574,4.02694776337728,100,982.7908020019531,5 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,1,374285.5654676295,6.041712643937134,62,3104.1543350219727,4 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,4,378759.2477216333,6.0447506876567285,56,2764.367919921875,4 +LN_COBYLA,Ipopt,0.0,1,368770.77971220895,5.154334896050398,100,1043.3742561340332,5 +LN_COBYLA,Ipopt,0.0,2,371232.7086212399,5.234096040616625,100,1079.748556137085,5 +LN_COBYLA,Ipopt,0.0,3,340062.9032832226,4.925229352545818,100,1151.005206823349,5 +LN_COBYLA,Ipopt,0.0,4,368212.45521732606,5.284610897966902,100,1045.847148180008,5 +LN_COBYLA,Ipopt,0.0,5,365393.8200764219,5.200453644061956,100,1072.3246850967407,5 +LN_COBYLA,Ipopt,0.0,6,336817.12414661905,4.700154258913552,100,1098.9553389549255,5 +LN_COBYLA,Ipopt,0.0,7,362356.2631309275,4.952592988627873,100,1035.2730069160461,5 +LN_COBYLA,Ipopt,0.0,8,364983.8689721409,5.103900633505085,100,1041.5780990123749,5 +LN_COBYLA,Ipopt,0.0,9,345868.855240273,4.772237635224206,100,1089.8890671730042,5 +LN_COBYLA,Ipopt,0.0,10,375547.58596656244,5.284457465353491,100,1023.4139959812164,5 +LN_NELDERMEAD,Ipopt,0.0,1,368770.7797122093,5.154334896050398,100,1029.9820940494537,5 +LN_NELDERMEAD,Ipopt,0.0,2,371232.7086212405,5.234096040616625,100,1155.5187888145447,5 +LN_NELDERMEAD,Ipopt,0.0,3,340062.90328358323,4.931194237636005,100,1123.0712430477142,5 +LN_NELDERMEAD,Ipopt,0.0,4,368212.4552173274,5.284610897966902,100,1079.9773318767548,5 +LN_NELDERMEAD,Ipopt,0.0,5,365393.82007642306,5.200453644061957,100,1069.9418380260468,5 +LN_NELDERMEAD,Ipopt,0.0,6,336817.12414661655,4.700154258913552,100,1003.5985488891602,5 +LN_NELDERMEAD,Ipopt,0.0,7,362356.26313091733,4.952592988627873,100,1048.8882961273193,5 +LN_NELDERMEAD,Ipopt,0.0,8,364983.86897214066,5.103900633505085,100,969.1135029792786,5 +LN_NELDERMEAD,Ipopt,0.0,9,345868.8552402731,4.772237635224206,100,1127.8717539310455,5 +LN_NELDERMEAD,Ipopt,0.0,10,375547.5859665647,5.2844574653534915,100,1156.0585510730743,5 +LD_CCSAQ,Ipopt,nothing,1,374287.53230167157,6.041903274738079,61,3224.8785610198975,4 +LD_CCSAQ,Ipopt,nothing,2,386275.29943069845,5.983625043111539,60,2969.909969806671,4 +LD_CCSAQ,Ipopt,nothing,3,357828.6772629269,6.32449112674849,61,2903.837574005127,4 +LD_CCSAQ,Ipopt,nothing,4,378759.3506686415,6.044780871279457,59,2952.851217985153,4 +LD_CCSAQ,Ipopt,nothing,5,377585.18819117936,6.187277191005196,61,3007.517009973526,4 +LD_CCSAQ,Ipopt,nothing,6,357932.69040892785,5.856561795392329,61,3102.083848953247,4 +LD_CCSAQ,Ipopt,nothing,7,365780.49135180935,5.613264808250105,59,3015.02862405777,4 +LD_CCSAQ,Ipopt,nothing,8,345286.2376051097,5.821766834284973,58,2842.756448984146,4 +LD_CCSAQ,Ipopt,nothing,9,352682.0504250703,6.401288787080837,60,3061.2380788326263,4 +LD_CCSAQ,Ipopt,nothing,10,364185.8161657268,6.141413958902178,59,2888.3543219566345,4 +LN_NEWUOA_BOUND,Ipopt,0.0,1,296056.5240687083,3.865897697776273,100,1021.5894551277161,5 +LN_NEWUOA_BOUND,Ipopt,0.0,2,297281.23558716115,3.9125716510869886,100,1070.2945890426636,5 +LN_NEWUOA_BOUND,Ipopt,0.0,3,275479.41710101644,3.7330186146693944,100,1144.276999950409,5 +LN_NEWUOA_BOUND,Ipopt,0.0,4,305627.4761332325,4.030881678722003,100,1118.0375318527222,5 +LN_NEWUOA_BOUND,Ipopt,0.0,5,293247.731179087,3.8324415568176824,100,1148.7592511177063,5 +LN_NEWUOA_BOUND,Ipopt,0.0,6,270883.3068218491,4.22518758895544,100,986.9172720909119,5 +LN_NEWUOA_BOUND,Ipopt,0.0,7,293612.63005215995,5.289572893298926,100,1073.409821987152,5 +LN_NEWUOA_BOUND,Ipopt,0.0,8,299671.27750323794,3.876275233464007,100,1015.0241770744324,5 +LN_NEWUOA_BOUND,Ipopt,0.0,9,284526.2713338175,6.1455362066452155,100,1145.198306798935,5 +LN_NEWUOA_BOUND,Ipopt,0.0,10,304363.3091196187,3.998324157548252,100,1121.8096508979797,5 +LD_SLSQP,Ipopt,nothing,1,367656.99680260976,6.309903805593335,100,4979.890875816345,5 +LD_SLSQP,Ipopt,nothing,2,297264.6518479442,3.9125716510869886,1,64.51116108894348,4 +LD_SLSQP,Ipopt,nothing,3,275474.4664725441,3.7330186146693944,1,60.47740888595581,4 +LD_SLSQP,Ipopt,nothing,4,305611.4541170079,4.030881678722003,1,64.25777292251587,4 +LD_SLSQP,Ipopt,nothing,5,377031.3091568417,6.279839303967756,100,5056.611438989639,5 +LD_SLSQP,Ipopt,nothing,6,356265.1239745704,6.038739548828272,85,4219.168741941452,4 +LD_SLSQP,Ipopt,nothing,7,293608.42417688074,3.759529176498192,1,62.80996298789978,4 +LD_SLSQP,Ipopt,nothing,8,353997.8791784004,5.186601584731852,100,4799.876461029053,5 +LD_SLSQP,Ipopt,nothing,9,284516.79877157696,3.779076104265464,1,62.72271800041199,4 +LD_SLSQP,Ipopt,nothing,10,363296.23298497265,6.016711778607968,100,5127.965346813202,5 diff --git a/results/strategic_bidding_nlopt_pglib_opf_case2000_goc.pdf b/results/strategic_bidding_nlopt_pglib_opf_case2000_goc.pdf new file mode 100644 index 00000000..1c9b6b18 Binary files /dev/null and b/results/strategic_bidding_nlopt_pglib_opf_case2000_goc.pdf differ diff --git a/results/strategic_bidding_nlopt_pglib_opf_case2868_rte.m.csv b/results/strategic_bidding_nlopt_pglib_opf_case2868_rte.m.csv new file mode 100644 index 00000000..e4bedfa6 --- /dev/null +++ b/results/strategic_bidding_nlopt_pglib_opf_case2868_rte.m.csv @@ -0,0 +1,74 @@ +solver_upper,solver_lower,Δp,seed,profit,market_share,num_evals,time,status +LD_MMA,Ipopt,nothing,1,629583.9288212913,6.890645774703093,58,3505.258134126663,4 +LD_MMA,Ipopt,nothing,2,627113.0582729869,6.977007604840815,57,3446.124926805496,4 +LD_MMA,Ipopt,nothing,3,503483.0698924671,9.437976029449642,31,1844.9283549785614,4 +LD_MMA,Ipopt,nothing,4,660150.8766179021,6.745448274791942,54,3128.340518951416,4 +LD_MMA,Ipopt,nothing,5,515708.90176392137,8.192006367082637,55,3282.953973054886,4 +LD_MMA,Ipopt,nothing,6,593330.5314686802,6.607485304488751,51,3126.637830018997,4 +LD_MMA,Ipopt,nothing,7,641960.0279960671,6.677319587937669,56,3369.3064410686493,4 +LD_MMA,Ipopt,nothing,8,441400.4075766634,9.50172834519057,52,3100.031136035919,4 +LD_MMA,Ipopt,nothing,9,505161.25023751793,9.46362463803075,54,3583.3756749629974,4 +LD_MMA,Ipopt,nothing,10,493257.18844268494,9.283018737480724,47,2828.7756350040436,4 +LN_BOBYQA,Ipopt,0.0,1,398781.7018975408,2.681879516265251,100,1111.6657028198242,5 +LN_BOBYQA,Ipopt,0.0,2,424994.5800923622,2.8237977709425417,100,1057.730427980423,5 +LN_BOBYQA,Ipopt,0.0,3,413469.5047721841,2.7924051654318505,100,1075.452929019928,5 +LN_BOBYQA,Ipopt,0.0,4,408094.1512728784,2.703931462513121,100,1070.8178579807281,5 +LN_BOBYQA,Ipopt,0.0,5,391744.736736471,2.6342395759849984,100,1036.4698250293732,5 +LN_BOBYQA,Ipopt,0.0,6,405613.8253293786,2.747279007315161,100,1144.034110069275,5 +LN_BOBYQA,Ipopt,0.0,7,405356.3161075126,2.714075307402132,100,1039.9613230228424,5 +LN_BOBYQA,Ipopt,0.0,8,416907.8601131197,2.8107092203326007,100,1105.951191186905,5 +LN_BOBYQA,Ipopt,0.0,9,407432.24270854925,2.7734348694199,100,1056.5739369392395,5 +LN_BOBYQA,Ipopt,0.0,10,388927.06299417454,2.6444588672801075,100,1018.8784649372101,5 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,4,644985.796929955,8.088733324695841,38,2319.878350019455,4 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,5,604839.8342306449,6.224460851907855,65,3976.07225894928,4 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,8,579232.5206609542,5.863372546358108,56,3459.778429031372,4 +LN_COBYLA,Ipopt,0.0,1,473527.52936282655,3.2914449436356303,100,1078.0652148723602,5 +LN_COBYLA,Ipopt,0.0,2,453459.1354854593,3.0572512242586547,100,1023.7090928554535,5 +LN_COBYLA,Ipopt,0.0,3,489848.76781296195,3.4914091683973116,100,1097.9070301055908,5 +LN_COBYLA,Ipopt,0.0,4,483885.64886411006,3.3273901264543326,100,962.7221291065216,5 +LN_COBYLA,Ipopt,0.0,5,456512.7012026482,3.2580187147012483,100,1095.7954370975494,5 +LN_COBYLA,Ipopt,0.0,6,475336.5075895495,3.408706698951814,100,1169.6885168552399,5 +LN_COBYLA,Ipopt,0.0,7,484500.91024239425,3.413132685829475,100,989.1927120685577,5 +LN_COBYLA,Ipopt,0.0,8,486279.8287879631,3.457433477227418,100,1181.4076688289642,5 +LN_COBYLA,Ipopt,0.0,9,460900.76954934996,3.2901854106692205,100,1078.2053899765015,5 +LN_COBYLA,Ipopt,0.0,10,470602.82532957214,3.334392838176782,100,1000.6699910163879,5 +LN_NELDERMEAD,Ipopt,0.0,1,473527.5293628363,3.2914449436356303,100,1118.1400711536407,5 +LN_NELDERMEAD,Ipopt,0.0,2,453459.1354854287,3.0572512242586547,100,1039.071209192276,5 +LN_NELDERMEAD,Ipopt,0.0,3,489848.76781296177,3.4914091683973116,100,1075.2348010540009,5 +LN_NELDERMEAD,Ipopt,0.0,4,483885.6488640899,3.3273901264543326,100,1033.4427018165588,5 +LN_NELDERMEAD,Ipopt,0.0,5,456512.6996559689,3.2580187147012483,100,984.0106539726257,5 +LN_NELDERMEAD,Ipopt,0.0,6,475336.507588969,3.408706698951814,100,1242.2352900505066,5 +LN_NELDERMEAD,Ipopt,0.0,7,484500.9102423927,3.413132685829475,100,1054.6953508853912,5 +LN_NELDERMEAD,Ipopt,0.0,8,486279.82878796756,3.457433477227418,100,1140.80641913414,5 +LN_NELDERMEAD,Ipopt,0.0,9,460900.7695494116,3.2901854106692205,100,1067.6936299800873,5 +LN_NELDERMEAD,Ipopt,0.0,10,470602.82532957214,3.334392838176782,100,1035.6403188705444,5 +LD_CCSAQ,Ipopt,nothing,1,624324.240321726,7.229484587659512,59,3568.3671629428864,4 +LD_CCSAQ,Ipopt,nothing,2,633243.8458275096,7.202017625481673,52,3150.5890390872955,4 +LD_CCSAQ,Ipopt,nothing,3,604991.2995579528,6.851712961741274,56,3349.7063908576965,4 +LD_CCSAQ,Ipopt,nothing,4,654112.6526452895,6.849097981255742,55,3309.5860052108765,4 +LD_CCSAQ,Ipopt,nothing,5,598384.9908440018,6.875877332376444,47,2825.262242794037,4 +LD_CCSAQ,Ipopt,nothing,6,599138.5529587846,6.570572979767129,43,2549.805088043213,4 +LD_CCSAQ,Ipopt,nothing,7,633182.7971276582,6.8977303515942365,57,3401.1067910194397,4 +LD_CCSAQ,Ipopt,nothing,8,563108.9678337795,7.421622914249337,59,3577.2335109710693,4 +LD_CCSAQ,Ipopt,nothing,9,605788.7008906608,6.966615464773943,57,3441.102548122406,4 +LD_CCSAQ,Ipopt,nothing,10,560763.6032972333,6.715269444455808,59,3569.4729290008545,4 +LN_NEWUOA_BOUND,Ipopt,0.0,1,396694.5937481017,2.668257454596464,100,1160.3855929374695,5 +LN_NEWUOA_BOUND,Ipopt,0.0,2,423004.55900321295,112.7964266242699,100,1102.4114229679108,5 +LN_NEWUOA_BOUND,Ipopt,0.0,3,411612.76630639104,3.1219005757152187,100,1098.7946500778198,5 +LN_NEWUOA_BOUND,Ipopt,0.0,4,406143.30047178635,2.6898943004534064,100,975.0623390674591,5 +LN_NEWUOA_BOUND,Ipopt,0.0,5,389837.04905492143,2.621566390611279,100,1086.1601371765137,5 +LN_NEWUOA_BOUND,Ipopt,0.0,6,403540.38210378896,2.5112192725532276,100,1164.4824831485748,5 +LN_NEWUOA_BOUND,Ipopt,0.0,7,403425.04139754584,2.7011857301054225,100,1085.67134308815,5 +LN_NEWUOA_BOUND,Ipopt,0.0,8,415381.09212414495,2.7980494975362054,100,1140.2639210224152,5 +LN_NEWUOA_BOUND,Ipopt,0.0,9,404450.27723066194,2.762864452960727,100,1073.3985180854797,5 +LN_NEWUOA_BOUND,Ipopt,0.0,10,386770.63151524833,2.630424732075897,100,1050.3003928661346,5 +LD_SLSQP,Ipopt,nothing,1,396692.3285684554,2.668257454596464,1,79.38984894752502,4 +LD_SLSQP,Ipopt,nothing,2,422994.42245616694,2.8107528763948797,1,76.5017671585083,4 +LD_SLSQP,Ipopt,nothing,3,595256.5472788628,5.229886902134434,100,6124.816920995712,5 +LD_SLSQP,Ipopt,nothing,4,406133.88114358974,2.6898943004534064,1,75.31735301017761,4 +LD_SLSQP,Ipopt,nothing,5,599438.027210647,6.302735019143101,100,6119.690332889557,5 +LD_SLSQP,Ipopt,nothing,6,403529.5235328151,2.7339725059105686,1,77.51345705986023,4 +LD_SLSQP,Ipopt,nothing,7,635602.1987549539,6.280930279767297,100,6045.816938877106,5 +LD_SLSQP,Ipopt,nothing,8,573160.8856220098,5.296483320893043,100,6050.58663392067,5 +LD_SLSQP,Ipopt,nothing,9,602261.8160636103,6.863588314096251,100,6154.7616510391235,5 +LD_SLSQP,Ipopt,nothing,10,386764.4662304388,2.630424732075897,1,76.16425085067749,4 diff --git a/results/strategic_bidding_nlopt_pglib_opf_case2868_rte.m.pdf b/results/strategic_bidding_nlopt_pglib_opf_case2868_rte.m.pdf new file mode 100644 index 00000000..547e1973 Binary files /dev/null and b/results/strategic_bidding_nlopt_pglib_opf_case2868_rte.m.pdf differ diff --git a/results/strategic_bidding_nlopt_pglib_opf_case2869_pegase.csv b/results/strategic_bidding_nlopt_pglib_opf_case2869_pegase.csv new file mode 100644 index 00000000..d6eea0a1 --- /dev/null +++ b/results/strategic_bidding_nlopt_pglib_opf_case2869_pegase.csv @@ -0,0 +1,66 @@ +solver_upper,solver_lower,Δp,seed,profit,market_share,num_evals,time,status +LD_MMA,Ipopt,nothing,1,763095.0244814962,5.473044706849988,58,5223.248332023621,4 +LD_MMA,Ipopt,nothing,2,647615.6240026074,5.510982189555721,51,4742.689513921738,4 +LD_MMA,Ipopt,nothing,3,785222.7713789152,5.496784911947918,52,4898.633311986923,4 +LD_MMA,Ipopt,nothing,4,696144.0876179516,5.693404697323059,48,5097.842547178268,4 +LD_MMA,Ipopt,nothing,5,777237.1560935642,5.465085022928578,31,3426.131083011627,4 +LD_MMA,Ipopt,nothing,7,709407.5791562888,5.551524568785184,51,4731.625789880753,4 +LD_MMA,Ipopt,nothing,8,829592.6036594185,5.503763274506792,54,5420.210821151733,4 +LD_MMA,Ipopt,nothing,9,740503.4336833901,5.494251454803217,53,4863.309112071991,4 +LN_BOBYQA,Ipopt,0.0,2,385557.7831737955,1.0457149046685181,100,2048.12282204628,5 +LN_BOBYQA,Ipopt,0.0,3,384922.0106639474,1.0338236542586012,100,1741.6074459552765,5 +LN_BOBYQA,Ipopt,0.0,4,372987.8115548726,1.000955716353836,100,1879.1780288219452,5 +LN_BOBYQA,Ipopt,0.0,8,385599.3639342813,1.040553129624309,100,1853.8286619186401,5 +LN_BOBYQA,Ipopt,0.0,9,377596.75082265434,1.0280291628739298,100,1888.7991490364075,5 +LN_BOBYQA,Ipopt,0.0,10,364604.2245929031,0.9789858879518353,100,1799.6835939884186,5 +LN_COBYLA,Ipopt,0.0,2,455690.4622564906,1.2768479925360656,100,2056.48614692688,5 +LN_COBYLA,Ipopt,0.0,3,453681.9058440277,1.268292117138062,100,1788.0797748565674,5 +LN_COBYLA,Ipopt,0.0,4,445688.2818385775,1.243156638955702,100,1846.7739961147308,5 +LN_COBYLA,Ipopt,0.0,8,402582.56670932414,1.1039157409656242,100,1886.1198689937592,5 +LN_COBYLA,Ipopt,0.0,9,460929.08241540554,1.2936620440698983,100,1753.9974720478058,5 +LN_COBYLA,Ipopt,0.0,10,441021.88283406827,1.234313197882759,100,1746.5269379615784,5 +LN_COBYLA,Ipopt,0.0,1,445078.1423079036,1.2348499223505311,100,1642.233705997467,5 +LD_SLSQP,Ipopt,nothing,5,361092.8467074701,0.9705346029574599,1,328.6287798881531,4 +LD_CCSAQ,Ipopt,nothing,1,865777.6409374254,4.924741812075953,58,4722.888010025024,4 +LN_BOBYQA,Ipopt,0.0,6,377512.7432059537,1.0171404460731441,100,2210.409973859787,5 +LN_BOBYQA,Ipopt,0.0,7,372795.4990087037,1.0051094144550203,100,4360.791481018066,5 +LD_MMA,Ipopt,nothing,6,655217.2314237518,5.768599646056398,51,4790.7436101436615,4 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,4,946819.1086822012,4.2721327415139125,61,6084.313436985016,4 +LD_MMA,Ipopt,nothing,10,740397.4707648519,5.437243652642576,51,4644.582688808441,4 +LD_CCSAQ,Ipopt,nothing,2,958960.7569564288,4.9773509468261565,53,4771.052474975586,4 +LD_CCSAQ,Ipopt,nothing,3,881844.446174467,5.030087253722467,52,4612.008496999741,4 +LN_BOBYQA,Ipopt,0.0,1,370868.89833338564,0.99269379010979,100,1574.4580879211426,5 +LD_CCSAQ,Ipopt,nothing,4,852544.2212812195,4.937471040422679,57,6256.860535144806,4 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,8,886900.2692065428,4.522389761538202,85,9017.482391119003,4 +LD_CCSAQ,Ipopt,nothing,5,876558.4320741579,4.9579754656702315,60,6208.630883932114,4 +LD_CCSAQ,Ipopt,nothing,6,882604.859439062,4.949006378588906,57,5059.946820020676,4 +LD_CCSAQ,Ipopt,nothing,7,865919.200103611,4.936893793045543,56,4742.679541110992,4 +LD_SLSQP,Ipopt,nothing,4,939773.5724978196,3.712319507250259,100,9312.706056833267,5 +LD_SLSQP,Ipopt,nothing,7,370846.8680343175,1.000010615588849,1,119.07087993621826,4 +LD_SLSQP,Ipopt,nothing,2,981751.7547293808,3.976173171005669,100,9234.672760009766,5 +LD_SLSQP,Ipopt,nothing,8,880621.0068699884,3.2708947211254404,100,8346.3530189991,5 +LD_SLSQP,Ipopt,nothing,9,962220.4989043615,4.253120976146774,100,8124.440047979355,5 +LD_SLSQP,Ipopt,nothing,10,915427.2589328055,4.072749203205618,100,8140.769929885864,5 +LD_SLSQP,Ipopt,nothing,3,382965.20115181385,1.0285899909371885,1,109.31101393699646,4 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,3,894864.6797276961,4.501147755337157,44,3928.2089171409607,4 +LD_CCSAQ,Ipopt,nothing,8,843574.6079864681,4.987640110007977,47,3896.163267850876,4 +LD_CCSAQ,Ipopt,nothing,9,893301.5427800094,4.964824798862503,53,4455.742611169815,4 +LD_CCSAQ,Ipopt,nothing,10,870338.9921688911,4.91134686506627,60,4816.567409992218,4 +LD_SLSQP,Ipopt,nothing,1,963555.3278867417,4.030547695525361,100,8082.192088127136,5 +LN_NELDERMEAD,Ipopt,0.0,1,445078.1430213009,1.2348499223505311,100,1329.5679678916931,5 +LN_NELDERMEAD,Ipopt,0.0,2,455690.4622564971,1.2768479925360656,100,1422.1589450836182,5 +LN_NELDERMEAD,Ipopt,0.0,3,453681.89531589224,1.268292117138062,100,1389.1774249076843,5 +LN_NELDERMEAD,Ipopt,0.0,4,445688.28183857014,1.2431566389557018,100,1524.381493806839,5 +LN_NELDERMEAD,Ipopt,0.0,5,438022.6043508742,1.2197590396224085,100,1577.4441611766815,5 +LN_NELDERMEAD,Ipopt,0.0,7,452235.4254595,1.2748403654661529,100,2957.7306208610535,5 +LN_NELDERMEAD,Ipopt,0.0,8,402582.56670931267,1.1039157409656242,100,1325.344002008438,5 +LN_NELDERMEAD,Ipopt,0.0,9,460929.0832423285,1.2936620440698983,100,1377.7193779945374,5 +LN_NELDERMEAD,Ipopt,0.0,10,441021.8828340604,1.234313197882759,100,1383.0285909175873,5 +LN_NEWUOA_BOUND,Ipopt,0.0,1,369089.77692523686,0.9878201820710817,100,1370.1421120166779,5 +LN_NEWUOA_BOUND,Ipopt,0.0,2,383773.56707333704,1.0405736572886726,100,1459.6795320510864,5 +LN_NEWUOA_BOUND,Ipopt,0.0,3,382966.3199129843,1.0285899909371885,100,1412.3048560619354,5 +LN_NEWUOA_BOUND,Ipopt,0.0,4,371114.4143419697,0.9958303959944159,100,1539.8862872123718,5 +LN_NEWUOA_BOUND,Ipopt,0.0,7,370847.88196845306,1.000010615588849,100,1651.1666140556335,5 +LN_NEWUOA_BOUND,Ipopt,0.0,8,384438.4956870938,1.0358707175496766,100,1306.5130269527435,5 +LN_NEWUOA_BOUND,Ipopt,0.0,9,375783.6996468839,1.0228448016738099,100,1501.0024898052216,5 +LN_NEWUOA_BOUND,Ipopt,0.0,10,362687.60243511904,0.9738140647887583,100,1411.0102169513702,5 diff --git a/results/strategic_bidding_nlopt_pglib_opf_case2869_pegase.pdf b/results/strategic_bidding_nlopt_pglib_opf_case2869_pegase.pdf new file mode 100644 index 00000000..9b9f5455 Binary files /dev/null and b/results/strategic_bidding_nlopt_pglib_opf_case2869_pegase.pdf differ diff --git a/results/strategic_bidding_nlopt_pglib_opf_case300_ieee.csv b/results/strategic_bidding_nlopt_pglib_opf_case300_ieee.csv new file mode 100644 index 00000000..978dbdf9 --- /dev/null +++ b/results/strategic_bidding_nlopt_pglib_opf_case300_ieee.csv @@ -0,0 +1,101 @@ +solver_upper,solver_lower,Δp,seed,profit,market_share,num_evals,time,status +LD_MMA,Ipopt,nothing,1,96110.71882678397,51.409912646352794,55,112.35776281356812,1 +LD_MMA,Ipopt,nothing,2,169611.93004530968,7.801358484116504,50,69.22345805168152,1 +LD_MMA,Ipopt,nothing,3,178444.28138203637,15.922922168322174,49,65.88223600387573,1 +LD_MMA,Ipopt,nothing,4,194336.04203766666,12.075242297979434,44,56.839174032211304,1 +LD_MMA,Ipopt,nothing,5,176613.2149870944,12.422180711866769,62,85.81530904769897,1 +LD_MMA,Ipopt,nothing,6,65458.26361763022,26.589452398348463,35,50.73106598854065,1 +LD_MMA,Ipopt,nothing,7,79498.15073959171,26.38580354505366,42,58.235145807266235,1 +LD_MMA,Ipopt,nothing,8,74496.46392915075,24.967103266479842,39,51.71993613243103,1 +LD_MMA,Ipopt,nothing,9,181934.25932745318,7.077167283620686,68,96.30612301826477,1 +LD_MMA,Ipopt,nothing,10,150551.556318804,15.517238331302416,44,61.61980485916138,1 +LN_BOBYQA,Ipopt,0.0,1,244232.89513945184,10.81443490086525,100,72.240394115448,1 +LN_BOBYQA,Ipopt,0.0,2,188387.21949619512,6.913013068576153,100,115.08367490768433,1 +LN_BOBYQA,Ipopt,0.0,3,193988.9162058156,7.025273557559125,100,65.74386215209961,1 +LN_BOBYQA,Ipopt,0.0,4,150513.5414421936,4.749046588752977,100,68.16472101211548,1 +LN_BOBYQA,Ipopt,0.0,5,237987.9450137697,8.291928779369465,100,70.9854621887207,1 +LN_BOBYQA,Ipopt,0.0,6,232213.14261507624,7.731145691394789,100,76.82888698577881,1 +LN_BOBYQA,Ipopt,0.0,7,237795.28401277255,10.070900059768746,100,73.00460886955261,1 +LN_BOBYQA,Ipopt,0.0,8,224848.85898168833,7.89813901386081,100,66.71672105789185,1 +LN_BOBYQA,Ipopt,0.0,9,274612.1997313932,11.68412522749941,100,67.7688729763031,1 +LN_BOBYQA,Ipopt,0.0,10,173580.3273015472,8.560154430382207,100,71.89584803581238,1 +LD_LBFGS,Ipopt,nothing,1,221645.11888298934,11.401696814027988,26,39.10578203201294,1 +LD_LBFGS,Ipopt,nothing,2,141262.71448922117,26.311399334821562,26,38.36118006706238,1 +LD_LBFGS,Ipopt,nothing,3,87363.0693632695,24.253863875152994,69,95.3198938369751,1 +LD_LBFGS,Ipopt,nothing,4,223731.48169915998,15.140592494373312,28,38.29574799537659,1 +LD_LBFGS,Ipopt,nothing,5,160116.8615679772,22.860132440024824,24,33.951844930648804,1 +LD_LBFGS,Ipopt,nothing,6,92192.29686791127,22.910628369499896,23,34.021559953689575,1 +LD_LBFGS,Ipopt,nothing,7,235055.40643308955,13.34543966105306,55,79.17287993431091,1 +LD_LBFGS,Ipopt,nothing,8,108555.6621307289,23.11464351812243,33,44.80105805397034,1 +LD_LBFGS,Ipopt,nothing,9,160994.8478759525,24.081755403931364,39,58.67839002609253,1 +LD_LBFGS,Ipopt,nothing,10,190046.33552746775,10.594052014825742,60,84.64722919464111,1 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,1,198035.28078277098,7.88719158718009,37,52.36943316459656,1 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,2,145472.39213218822,26.120541699500915,49,74.90275287628174,1 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,3,87359.72330924864,24.251828921465503,82,122.74614095687866,1 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,4,223801.14490554525,15.217660857646734,34,47.43060493469238,1 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,5,160116.8615679772,22.860132440024824,25,35.79092884063721,1 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,6,92192.29686791127,22.910628369499896,24,35.60436296463013,1 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,7,253347.71115981837,11.747118556625795,95,138.99833989143372,1 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,8,108543.39782339208,23.10815061751523,46,62.6436607837677,1 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,9,161039.6722133498,24.06518542269386,65,96.75429606437683,1 +LD_TNEWTON_PRECOND_RESTART,Ipopt,nothing,10,190069.3908097993,10.536935293217182,75,104.32000207901001,1 +LN_COBYLA,Ipopt,0.0,1,176010.19169836925,5.677305066259365,100,72.03100514411926,1 +LN_COBYLA,Ipopt,0.0,2,127374.33544940113,3.8448240810590555,100,64.18761801719666,1 +LN_COBYLA,Ipopt,0.0,3,176150.12057850952,5.482315426210346,100,65.934485912323,1 +LN_COBYLA,Ipopt,0.0,4,91727.20714677343,2.8158531169098886,100,68.77439904212952,1 +LN_COBYLA,Ipopt,0.0,5,198362.8753737887,6.237497498689416,100,72.6327908039093,1 +LN_COBYLA,Ipopt,0.0,6,102785.83811041675,2.9832594145233964,100,70.44147896766663,1 +LN_COBYLA,Ipopt,0.0,7,187911.54789161313,5.893495185002765,100,73.67294597625732,1 +LN_COBYLA,Ipopt,0.0,8,148226.30580434453,4.7518733303087295,100,66.37015104293823,1 +LN_COBYLA,Ipopt,0.0,9,106883.42897620534,3.1732717624926936,100,68.65271401405334,1 +LN_COBYLA,Ipopt,0.0,10,105033.57597358119,3.1867592685601807,100,68.18176698684692,1 +LD_CCSAQ,Ipopt,nothing,1,238249.22192024256,13.368774486955296,53,77.45764899253845,1 +LD_CCSAQ,Ipopt,nothing,2,152988.0994922613,19.838917943421166,51,73.36208820343018,1 +LD_CCSAQ,Ipopt,nothing,3,54459.19634512048,32.066785148062614,53,79.32646608352661,1 +LD_CCSAQ,Ipopt,nothing,4,223330.3099661669,14.979136317486603,59,77.34801506996155,1 +LD_CCSAQ,Ipopt,nothing,5,162860.81792287217,22.377943402548755,54,74.73434710502625,1 +LD_CCSAQ,Ipopt,nothing,6,238090.28459416158,16.0151547078734,56,77.43051195144653,1 +LD_CCSAQ,Ipopt,nothing,7,154809.254065423,16.97647571581712,55,78.2419638633728,1 +LD_CCSAQ,Ipopt,nothing,8,186665.13050085382,13.90320875043978,54,73.58215618133545,1 +LD_CCSAQ,Ipopt,nothing,9,103830.54778112366,29.50307082154892,46,69.35896801948547,1 +LD_CCSAQ,Ipopt,nothing,10,190069.37436692233,10.537467135211093,56,75.66284680366516,1 +LD_SLSQP,Ipopt,nothing,1,202122.28641528555,2.267153756014537,100,131.90146803855896,1 +LD_SLSQP,Ipopt,nothing,2,199331.8254757841,2.5237751272297952,100,138.61620903015137,1 +LD_SLSQP,Ipopt,nothing,3,223750.4938407642,2.2014218689308698,100,134.73497080802917,1 +LD_SLSQP,Ipopt,nothing,4,209293.7385906039,2.4208293676255077,100,129.0216929912567,1 +LD_SLSQP,Ipopt,nothing,5,203328.36023899372,1.4604130279370071,97,135.48368501663208,1 +LD_SLSQP,Ipopt,nothing,6,231355.38333033636,2.0340602661021707,100,140.8066110610962,1 +LD_SLSQP,Ipopt,nothing,7,242495.6487788387,4.432860622983777,100,136.01385188102722,1 +LD_SLSQP,Ipopt,nothing,8,187208.17397587682,1.6937783334110903,100,127.55766701698303,1 +LD_SLSQP,Ipopt,nothing,9,190462.07305198105,2.353980987448582,100,133.07019090652466,1 +LD_SLSQP,Ipopt,nothing,10,184874.85776357268,2.6584992062203487,100,135.63582706451416,1 +LN_NELDERMEAD,Ipopt,0.0,1,142797.72464495822,5.161901992702432,100,80.62988901138306,5 +LN_NELDERMEAD,Ipopt,0.0,2,162906.03031494704,6.7877009221550555,100,82.2264289855957,5 +LN_NELDERMEAD,Ipopt,0.0,3,187053.40210264924,6.456494149939711,100,76.6717300415039,5 +LN_NELDERMEAD,Ipopt,0.0,4,122654.05499742832,4.8276997818123775,100,71.19857501983643,5 +LN_NELDERMEAD,Ipopt,0.0,5,181280.52587454475,5.6929535395111825,100,79.03099584579468,5 +LN_NELDERMEAD,Ipopt,0.0,6,143977.16792034736,5.214501855637245,100,79.06425786018372,5 +LN_NELDERMEAD,Ipopt,0.0,7,191823.29265087267,6.515167628903469,100,85.2777590751648,5 +LN_NELDERMEAD,Ipopt,0.0,8,129787.3066354752,4.8447342062259375,100,74.79859685897827,5 +LN_NELDERMEAD,Ipopt,0.0,9,125691.36379888249,4.614137285862181,100,77.564857006073,5 +LN_NELDERMEAD,Ipopt,0.0,10,115360.68662683672,3.744635559065837,100,75.56841897964478,5 +LN_NEWUOA_BOUND,Ipopt,0.0,1,114609.03285998756,3.5386258929295185,100,77.19337701797485,5 +LN_NEWUOA_BOUND,Ipopt,0.0,2,126877.25981097147,4.278105995840755,100,80.73210501670837,5 +LN_NEWUOA_BOUND,Ipopt,0.0,3,179484.56136182434,5.774421048335941,100,77.0150978565216,5 +LN_NEWUOA_BOUND,Ipopt,0.0,4,69624.62873078427,2.2157777785271127,100,72.36222410202026,5 +LN_NEWUOA_BOUND,Ipopt,0.0,5,185082.98766657148,6.148943205505296,100,80.79146790504456,5 +LN_NEWUOA_BOUND,Ipopt,0.0,6,159079.53750403182,5.2054481202767695,100,80.37452101707458,5 +LN_NEWUOA_BOUND,Ipopt,0.0,7,244099.06401609167,10.392530675873035,100,86.66114091873169,5 +LN_NEWUOA_BOUND,Ipopt,0.0,8,41931.002594318015,1.18971351086284,100,75.69124221801758,5 +LN_NEWUOA_BOUND,Ipopt,0.0,9,186785.7114113925,6.29186177667848,100,78.50683999061584,5 +LN_NEWUOA_BOUND,Ipopt,0.0,10,30108.3054170366,0.9800887548124835,100,79.19660210609436,5 +G_MLSL_LDS,Ipopt,0.0,1,105853.85021551093,40.9637678014902,100,85.87584209442139,5 +G_MLSL_LDS,Ipopt,0.0,2,161395.64760316032,35.422872229648874,100,87.4148759841919,5 +G_MLSL_LDS,Ipopt,0.0,3,85866.78664958006,35.958121523824374,100,90.04950094223022,5 +G_MLSL_LDS,Ipopt,0.0,4,104546.12177641856,34.056326861683445,100,76.05526995658875,5 +G_MLSL_LDS,Ipopt,0.0,5,132577.46325762256,40.43476393507388,100,81.02201795578003,5 +G_MLSL_LDS,Ipopt,0.0,6,143977.16792034736,5.214501855637245,100,80.43115305900574,5 +G_MLSL_LDS,Ipopt,0.0,7,171313.9448243225,5.7162965290433005,100,83.96484804153442,5 +G_MLSL_LDS,Ipopt,0.0,8,129787.3066354752,4.8447342062259375,100,72.29825091362,5 +G_MLSL_LDS,Ipopt,0.0,9,125391.20812140916,35.508563827663444,100,92.51962304115295,5 +G_MLSL_LDS,Ipopt,0.0,10,145739.28083574225,35.126302700293266,100,83.76881694793701,5 diff --git a/results/strategic_bidding_nlopt_pglib_opf_case300_ieee.pdf b/results/strategic_bidding_nlopt_pglib_opf_case300_ieee.pdf new file mode 100644 index 00000000..78198a5f Binary files /dev/null and b/results/strategic_bidding_nlopt_pglib_opf_case300_ieee.pdf differ diff --git a/results/time_summary.csv b/results/time_summary.csv new file mode 100644 index 00000000..5fb5779f --- /dev/null +++ b/results/time_summary.csv @@ -0,0 +1,6 @@ +Best Improvement AVG,Best Improvement STD,Overall Improvement AVG,Overall Improvement STD,case +19.0,24.8327740429189,41.260000000000005,8.133360655691813,pglib_opf_case300_ieee +41.0,25.342103744997615,49.63166666666667,16.57324269588755,pglib_opf_case1354_pegase +23.7,22.7305472388541,42.50833333333333,16.63823732725389,pglib_opf_case2869_pegase +36.5,12.868998061663973,44.1,16.24431524496425,pglib_opf_case2000_goc +44.4,5.601587076693333,48.25,16.796714846351048,pglib_opf_case2868_rte.m diff --git a/slurm.jl b/slurm.jl new file mode 100644 index 00000000..980105f0 --- /dev/null +++ b/slurm.jl @@ -0,0 +1,19 @@ +# +using Pkg +Pkg.activate(@__DIR__) +try + + using Distributed, ClusterManagers +catch + Pkg.add("ClusterManagers") + Pkg.checkout("ClusterManagers") +end + +using Distributed, ClusterManagers + +np = 5 # +addprocs(SlurmManager(np), job_file_loc = ARGS[1]) + +println("We are all connected and ready.") + +include(ARGS[2]) \ No newline at end of file diff --git a/strategic_bidding.jl b/strategic_bidding.jl new file mode 100644 index 00000000..c566948b --- /dev/null +++ b/strategic_bidding.jl @@ -0,0 +1,270 @@ +using Distributed + +# Add worker processes if needed +# addprocs() + +@everywhere pkg_path = @__DIR__ + +@everywhere import Pkg + +@everywhere Pkg.activate(pkg_path) + +@everywhere Pkg.instantiate() + +@everywhere using JuMP, NLopt +@everywhere using DataFrames, CSV +@everywhere using SparseArrays +@everywhere using LinearAlgebra +@everywhere using Ipopt +@everywhere using Random +@everywhere using Logging + +@everywhere include("nlp_utilities.jl") +@everywhere include("nlp_utilities_test.jl") +@everywhere include("opf.jl") + +################################################ +# Strategic bidding test +################################################ + +# Parameters +@everywhere max_eval = 100 +@everywhere solver_lower_name = "Ipopt" +@everywhere casename = "pglib_opf_case300_ieee" # "pglib_opf_case300_ieee" "pglib_opf_case1354_pegase" "pglib_opf_case2869_pegase" "pglib_opf_case2000_goc" "pglib_opf_case2868_rte" +@everywhere save_file_name = "results/strategic_bidding_nlopt_$(casename)" +@everywhere save_file = save_file_name * ".csv" + +@everywhere data = make_basic_network(pglib(casename)) + +# #### test Range Evaluation +# Random.seed!(1) +# data = build_bidding_opf_model(casename) +# # offer = rand(length(data["bidding_lmps"])) +# # offer = zeros(length(data["bidding_lmps"])); offer[1] = 1.0; offer[2] = 1.0; offer[3] = 1.0; offer[4] = 1.0 +# offer = ones(length(data["bidding_lmps"])) +# profits = [] +# start_time = time() +# for val = 0.0:0.1:12 +# set_parameter_value.(data["bid"], val * offer) +# JuMP.optimize!(data["model"]) +# push!(profits, dot(-dual.(data["bidding_lmps"]), value.(data["bidding_generators_dispatch"]))) +# end +# end_time = time() +#### end test + +# #### test SB +# solver_upper = :LD_MMA # :LD_MMA :LN_BOBYQA +# Random.seed!(1234) +# start_time = time() +# profit, num_evals, trace, market_share, ret = test_bidding_nlopt(casename; percen_bidding_nodes=0.1, Δp=nothing, solver_lower=solver_lower, solver_upper=solver_upper, max_eval=2) +# end_time = time() +#### end test + +# Experiments +seeds = collect(1:10) + +experiements = Dict( + :LD_MMA => [nothing], + :LN_BOBYQA => [0.0], + :LD_CCSAQ => [nothing], + :LD_SLSQP => [nothing], + # :LD_LBFGS => [nothing], + :LD_TNEWTON_PRECOND_RESTART => [nothing], + :LN_COBYLA => [0.0], + :LN_NELDERMEAD => [0.0], + :LN_NEWUOA_BOUND => [0.0], + # :G_MLSL_LDS => [0.0], +) + +@everywhere res = Dict( + :FORCED_STOP => -5, + :ROUNDOFF_LIMITED => -4, + :OUT_OF_MEMORY => -3, + :INVALID_ARGS => -2, + :FAILURE => -1, + :SUCCESS => 1, + :STOPVAL_REACHED => 2, + :FTOL_REACHED => 3, + :XTOL_REACHED => 4, + :MAXEVAL_REACHED => 5, + :MAXTIME_REACHED => 6 +) + +# Prepare the list of experiments +_experiments = [ + (string(solver_upper), solver_lower_name, string(Δp), seed) + for (solver_upper, Δp_values) in experiements + for Δp in Δp_values + for seed in seeds +] + +# Check already executed experiments +if isfile(save_file) + old_results = CSV.read(save_file, DataFrame) + _experiments = setdiff(_experiments, [ + (string(row.solver_upper), string(row.solver_lower), string(row.Δp), row.seed) + for row in eachrow(old_results) + ]) +else + open(save_file, "w") do f + write(f, "solver_upper,solver_lower,Δp,seed,profit,market_share,num_evals,time,status\n") + end +end + +@everywhere function run_experiment(_solver_upper, _Δp, seed, id) + _data = deepcopy(data) + solver_upper = Symbol(_solver_upper) + Δp = _Δp == "nothing" ? nothing : parse(Float64, _Δp) + try + Random.seed!(seed) + start_time = time() + profit, num_evals, trace, market_share, ret = test_bidding_nlopt( + _data; percen_bidding_nodes=0.1, Δp=Δp, solver_upper=solver_upper, max_eval=max_eval + ) + end_time = time() + ret = res[ret] + if ret < 0 + @warn "Solver $(solver_upper) failed with seed $(seed)" + return nothing + else + df = DataFrame( + solver_upper = [string(solver_upper)], + solver_lower = [solver_lower_name], + Δp = [string(Δp)], + seed = [seed], + profit = [profit], + market_share = [market_share], + num_evals = [num_evals], + time = [end_time - start_time], + status = [ret] + ) + return df + end + catch e + @warn "Solver $(solver_upper) failed with seed $(seed)" e + return nothing + end +end + +# Run experiments using pmap +experiments_list = [ + (_solver_upper, _Δp, seed, id) + for (id, (_solver_upper, _, _Δp, seed)) in enumerate(_experiments) +] + +results = pmap( + experiment -> run_experiment(experiment...), + experiments_list +) + +# Filter out failed experiments +results = filter(!isnothing, results) + +# Combine results into a DataFrame +if !isempty(results) + all_results = vcat(results...) + # Append to existing results if any + if isfile(save_file) + old_results = CSV.read(save_file, DataFrame) + combined_results = vcat(old_results, all_results) + else + combined_results = all_results + end + # Save combined results + CSV.write(save_file, combined_results) +else + @info "No new results" +end + +# Proceed with plotting or further analysis as needed + + +# Plot results: One scatter point per solver_upper +# - Each solver_upper should be an interval with the mean and std of the results per seed +# - x-axis: market_share | y-axis: (profit / max_profit) * 100 +# - color: one per solver_upper +# - shape: one per Δp + +using Plots +using Statistics + +results = CSV.read(save_file, DataFrame) + +ignore_solver_upper = [:LD_LBFGS; :G_MLSL_LDS] + +results = results[[!(Symbol(sv_up) ∈ ignore_solver_upper) for sv_up in results.solver_upper], :] + +maximum_per_seed = [maximum(results.profit[results.seed .== seed]) for seed in seeds] + +results.gap = (maximum_per_seed[results.seed] - results.profit) * 100 ./ maximum_per_seed[results.seed] + +# use combined groupby +results_d = combine(groupby(results, [:solver_upper, :Δp]), :gap => mean, :market_share => mean, :gap => std, :num_evals => mean, :time => mean) + +# plot +markers = [] +for Δp in results_d.Δp + if Δp == "nothing" + push!(markers, :circle) + else + push!(markers, :diamond) + end +end +plt = scatter(results_d.market_share_mean, results_d.gap_mean, group=results_d.solver_upper, yerr=results_d.gap_std/2, + xlabel="Bid Volume (% of Market Share)", ylabel="Optimality Gap (%)", legend=:outertopright, marker=markers +) + +# save +savefig(plt, "results/strategic_bidding_nlopt_$(casename).pdf") + +# save a csv and plot for all cases the percentage gain (or loss) the best gradient-solver (Δp=="nothing") is when compared to the best non-gradient solver (Δp!=nothing) (inform mean and std) +using Plots +using Statistics +using DataFrames, CSV +seeds = collect(1:10) +cases = ["pglib_opf_case300_ieee", "pglib_opf_case1354_pegase", "pglib_opf_case2869_pegase", "pglib_opf_case2000_goc", "pglib_opf_case2868_rte.m"] +ignore_solver_upper = [:LD_LBFGS; :G_MLSL_LDS] +avg_gains = Array{Float64}(undef, length(cases)) +std_gains = Array{Float64}(undef, length(cases)) +avg_gains_all = Array{Float64}(undef, length(cases)) +std_gains_all = Array{Float64}(undef, length(cases)) +for (i, casename) in enumerate(cases) + save_file = "results/strategic_bidding_nlopt_$(casename).csv" + results = CSV.read(save_file, DataFrame) + results = results[[!(Symbol(sv_up) ∈ ignore_solver_upper) for sv_up in results.solver_upper], :] + results_d = combine(groupby(results, [:solver_upper, :Δp]), :profit => mean) + solver_upper_grad = results_d[results_d.Δp .== "nothing", :].solver_upper + solver_upper_nograd = results_d[results_d.Δp .!= "nothing", :].solver_upper + best_solver_upper_grad = solver_upper_grad[argmax(results_d[results_d.Δp .== "nothing", :].profit_mean)] + best_solver_upper_nograd = solver_upper_nograd[argmax(results_d[results_d.Δp .!= "nothing", :].profit_mean)] + mean_per_seed_non_gradient = [mean(results.profit[(results.seed .== seed) .& (results.Δp .!= "nothing")]) for seed in seeds] + mean_per_seed_gradient = [mean(results.profit[(results.seed .== seed) .& (results.Δp .== "nothing")]) for seed in seeds] + gains = [] + for seed in seeds + profit_grad = results[(results.seed .== seed) .* (results.solver_upper .== best_solver_upper_grad), :] + if isempty(profit_grad) + continue + end + profit_grad = profit_grad.profit[1] + profit_nograd = results[(results.seed .== seed) .* (results.solver_upper .== best_solver_upper_nograd), :] + if isempty(profit_nograd) + continue + end + profit_nograd = profit_nograd.profit[1] + push!(gains, (profit_grad - profit_nograd) * 100 / profit_nograd) + end + avg_gains[i] = mean(gains) + try + std_gains[i] = std(gains) + catch e + std_gains[i] = 0.0 + end + gains_all = (mean_per_seed_gradient .- mean_per_seed_non_gradient) * 100 ./ mean_per_seed_non_gradient + avg_gains_all[i] = mean(gains_all) + std_gains_all[i] = std(gains_all) +end + +# save +df = DataFrame(casename=cases, Best_Solver_Improvement_AVG=avg_gains, Best_Solver_Improvement_STD=std_gains, All_Solvers_Improvement_AVG=avg_gains_all, All_Solvers_Improvement_STD=std_gains_all) + + diff --git a/tables.jl b/tables.jl new file mode 100644 index 00000000..d709c7e2 --- /dev/null +++ b/tables.jl @@ -0,0 +1,66 @@ +using DataFrames, Statistics, CSV + +function compute_improvements(df::DataFrame) + # Get all unique seeds in the dataframe + seeds = unique(df.seed) + + best_improvements = Float64[] + overall_improvements = Float64[] + + for s in seeds + # Filter data for the current seed + df_seed = df[df.seed .== s, :] + + # Select rows for LD and LN methods. + # We assume that the 'solver_upper' column indicates "LD" or "LN" at the beginning. + ld = df_seed[startswith.(df_seed.solver_upper, "LD"), :] + ln = df_seed[startswith.(df_seed.solver_upper, "LN"), :] + + # Skip this seed if one of the groups is missing + if nrow(ld) == 0 || nrow(ln) == 0 + continue + end + + # For the best improvement: select the row with the best (highest) profit from each group. + ld_best = ld[findmax(ld.profit)[2], :] + ln_best = ln[findmax(ln.profit)[2], :] + + # Compute percentage improvement: how much lower LD's evaluations are compared to LN's. + best_impr = ((ln_best.num_evals - ld_best.num_evals) / ln_best.num_evals) * 100 + push!(best_improvements, best_impr) + + # For overall improvement: compute average num_evals for LD and LN. + ld_avg = mean(ld.num_evals) + ln_avg = mean(ln.num_evals) + overall_impr = ((ln_avg - ld_avg) / ln_avg) * 100 + push!(overall_improvements, overall_impr) + end + + # Compute overall averages and standard deviations across seeds. + best_avg = mean(best_improvements) + best_std = std(best_improvements) + overall_avg = mean(overall_improvements) + overall_std = std(overall_improvements) + + # Return the results as a one-row DataFrame. + return DataFrame( + "Best Improvement AVG" => best_avg, + "Best Improvement STD" => best_std, + "Overall Improvement AVG" => overall_avg, + "Overall Improvement STD" => overall_std, + ) +end + +# Example usage +prefix = "strategic_bidding_nlopt_" +case_names = ["pglib_opf_case300_ieee" "pglib_opf_case1354_pegase" "pglib_opf_case2869_pegase" "pglib_opf_case2000_goc" "pglib_opf_case2868_rte.m"] +df = DataFrame() +for case in case_names + file = "results/$(prefix)$(case).csv" + df_case = CSV.read(file, DataFrame) + df_results = compute_improvements(df_case) + df_results.case = [case] + append!(df, df_results) +end + +CSV.write("results/time_summary.csv", df) \ No newline at end of file diff --git a/test_jump.jl b/test_jump.jl new file mode 100644 index 00000000..abc4e974 --- /dev/null +++ b/test_jump.jl @@ -0,0 +1,175 @@ +using JuMP +using Ipopt + +model = JuMP.Model(Ipopt.Optimizer) + +@variable(model, x) +@variable(model, y) + +@objective(model, Min, x + y) + +con1_ = @constraint(model, x >= 1) + +con2_ = @constraint(model, -y <= -2) + +JuMP.optimize!(model) + +model.moi_backend.optimizer.model.inner.mult_g +model.moi_backend.optimizer.model.inner.g + +dual(con1_) + +dual(con2_) + +####### + +model = JuMP.Model(Ipopt.Optimizer) + +@variable(model, x) +@variable(model, y) +@variable(model, s1) +@variable(model, s2) + + +@objective(model, Min, x + y) + +con1 = @constraint(model, x - 1 - s1 == 0) + +con2 = @constraint(model, -y + 2 - s2 == 0) + +con1s = @constraint(model, s1 >= 0) +con2s = @constraint(model, s2 <= 0) + +JuMP.optimize!(model) + +dual(con1) +dual(con2) +dual(con1s) +dual(con2s) + +####### + +model = JuMP.Model(Ipopt.Optimizer) + +@variable(model, x) +@variable(model, y) + +@objective(model, Max, x + y) + +con1_ = @constraint(model, -x >= -1) + +con2_ = @constraint(model, y <= 2) + +JuMP.optimize!(model) + +model.moi_backend.optimizer.model.inner.mult_g +model.moi_backend.optimizer.model.inner.g + +dual(con1_) + +dual(con2_) + +####### + +model = JuMP.Model(Ipopt.Optimizer) + +@variable(model, 0 <= x <= 1) +@variable(model, 0 <= y <= 1) + +@objective(model, Max, x - y) + +JuMP.optimize!(model) + +model.moi_backend.optimizer.model.inner.mult_x_L +model.moi_backend.optimizer.model.inner.mult_x_U + +dual.(LowerBoundRef(y)) + +dual.(UpperBoundRef(x)) + + +####### + +model = JuMP.Model(Ipopt.Optimizer) + +@variable(model, 0 <= x <= 1) +@variable(model, 0 <= y <= 1) + +@objective(model, Min, x - y) + +JuMP.optimize!(model) + +model.moi_backend.optimizer.model.inner.mult_x_L + +model.moi_backend.optimizer.model.inner.mult_x_U + +dual.(LowerBoundRef(x)) + +dual.(UpperBoundRef(y)) + +####### + +model = JuMP.Model(Ipopt.Optimizer) + +@variable(model, x) +@variable(model, y <= 2) +@variable(model, z >= 3) +@variable(model, w) +@variable(model, u) + +@objective(model, Min, x - y + z + w - u) + +con1 = @constraint(model, x - 1 == 0) + +con1s = @constraint(model, w - 4 >= 0) +con2s = @constraint(model, u - 5 <= 0) + +JuMP.optimize!(model) + +JuMP.termination_status(model) + +model.moi_backend.optimizer.model.inner.mult_x_L +model.moi_backend.optimizer.model.inner.mult_x_U + +model.moi_backend.optimizer.model.inner.mult_g +model.moi_backend.optimizer.model.inner.g + +dual.(UpperBoundRef(y)) +dual.(LowerBoundRef(z)) +dual(con1) +dual(con1s) +dual(con2s) + + +####### Max + +model = JuMP.Model(Ipopt.Optimizer) + +@variable(model, x) +@variable(model, y <= 2) +@variable(model, z >= 3) +@variable(model, w) +@variable(model, u) + +@objective(model, Max, - x + y - z - w + u) + +con1 = @constraint(model, x - 1 == 0) +# s = w - 4 -> s >= 0 λ +con1s = @constraint(model, w - 4 >= 0) +con2s = @constraint(model, u - 5 <= 0) + +JuMP.optimize!(model) + +JuMP.termination_status(model) + +model.moi_backend.optimizer.model.inner.mult_x_L +model.moi_backend.optimizer.model.inner.mult_x_U + +model.moi_backend.optimizer.model.inner.mult_g +model.moi_backend.optimizer.model.inner.g + +dual.(UpperBoundRef(y)) +dual.(LowerBoundRef(z)) +dual(con1) +dual(con1s) +dual(con2s) \ No newline at end of file diff --git a/testing_barrier.jl b/testing_barrier.jl new file mode 100644 index 00000000..a6a7bcd8 --- /dev/null +++ b/testing_barrier.jl @@ -0,0 +1,78 @@ +################################################ +# sIPOPT test +################################################ + +using JuMP +using SparseArrays +using LinearAlgebra +using Ipopt +using Random +# using MadNLP +# using KNITRO + +include("nlp_utilities.jl") +include("nlp_utilities_test.jl") +include("opf.jl") + +# No Fix + +test_compute_optimal_hess_jacobian() + +test_compute_derivatives_Finite_Diff(DICT_PROBLEMS_no_cc, false) + +test_compute_derivatives_Analytical(DICT_PROBLEMS_Analytical_no_cc) + +# Fix and Relax + +test_compute_derivatives() + +test_compute_derivatives_Finite_Diff(DICT_PROBLEMS_cc, true) + +test_compute_derivatives_Analytical(DICT_PROBLEMS_Analytical_cc) + + +################################################ +# Moonlander test +################################################ + +include("moonlanding.jl") +using HSL_jll + +Isp_0 = 0.3 +model, tf, Isp = moonlander_JMP(_I=Isp_0) +set_optimizer(model, optimizer_with_attributes(Ipopt.Optimizer, + "print_level" => 0, + "hsllib" => HSL_jll.libhsl_path, + "linear_solver" => "MA57" +)) +JuMP.optimize!(model) +termination_status = JuMP.termination_status(model) +tf_0 = value(tf) + +function compute_sentitivity_for_var(model; _primal_vars=[], _cons=[], Δp=nothing) + primal_vars = all_primal_vars(model) + params = all_params(model) + vars_idx = [findall(x -> x == i, primal_vars)[1] for i in _primal_vars] + evaluator, cons = create_evaluator(model; x=[primal_vars; params]) + num_cons = length(cons) + cons_idx = [findall(x -> x == i, cons)[1] for i in _cons] + leq_locations, geq_locations = find_inequealities(cons) + num_ineq = length(leq_locations) + length(geq_locations) + num_primal = length(primal_vars) + + Δs, sp = compute_sensitivity(evaluator, cons, Δp; primal_vars=primal_vars, params=params) + return Δs[vars_idx, :], Δs[num_primal+num_ineq.+cons_idx, :], sp[vars_idx, :], sp[num_primal+num_ineq.+cons_idx, :] +end + +Δp=nothing + +Δs_primal, Δs_dual, sp_primal, sp_dual = compute_sentitivity_for_var(model; _primal_vars=[tf], Δp=Δp) +Δp = 0.001 +tf_p = tf_0 + Δs_primal[1] * Δp + +set_parameter_value(Isp, Isp_0+Δp) +JuMP.optimize!(model) +termination_status = JuMP.termination_status(model) +tf_p_true = value(tf) + +@test tf_p ≈ tf_p_true rtol=1e-3 \ No newline at end of file