From d1f1488cfc7ad23189b774628ee293715a93da77 Mon Sep 17 00:00:00 2001 From: thorek1 Date: Tue, 17 Oct 2023 08:48:36 +0200 Subject: [PATCH] got 1 period optim going --- Project.toml | 2 + src/MacroModelling.jl | 37 +++- src/macros.jl | 10 + src/structures.jl | 3 + test/amend_parser.jl | 476 ++++++++++++++++++++++++++++++++++++++++-- 5 files changed, 500 insertions(+), 28 deletions(-) diff --git a/Project.toml b/Project.toml index 802b37e4..dc1541e9 100644 --- a/Project.toml +++ b/Project.toml @@ -14,6 +14,7 @@ DynarePreprocessor_jll = "23afba7c-24e5-5ee2-bc2c-b42e07f0492a" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" ImplicitDifferentiation = "57b37032-215b-411a-8a7c-41a003a55207" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -30,6 +31,7 @@ RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" SpeedMapping = "f1835b91-879b-4a3f-a438-e4baacf14412" +StatusSwitchingQP = "73a7d86c-03e3-42fb-a0f5-70ef6cda08f4" Subscripts = "2b7f82d5-8785-4f63-971e-f18ddbeb808e" SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index a17d8b6f..9cd361cb 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -9,6 +9,7 @@ import SpecialFunctions: erfcinv, erfc import SymPyPythonCall as SPyPyC import Symbolics import ForwardDiff as β„± +import JuMP: AffExpr # import Zygote import SparseArrays: SparseMatrixCSC, SparseVector, AbstractSparseArray#, sparse, spzeros, droptol!, sparsevec, spdiagm, findnz#, sparse! import LinearAlgebra as β„’ @@ -124,6 +125,24 @@ Base.show(io::IO, 𝓂::β„³) = println(io, +function parse_obc_shock_bounds(expr::Expr) + # Determine the order of the shock and bound in the expression + shock_first = isa(expr.args[2], Symbol) + + # Extract the shock and bound from the expression + shock = shock_first ? expr.args[2] : expr.args[3] + bound_expr = shock_first ? expr.args[3] : expr.args[2] + + # Evaluate the bound expression to get a numerical value + bound = eval(bound_expr) |> Float64 + + # Determine whether the bound is a lower or upper bound based on the comparison operator and order + is_upper_bound = (expr.args[1] in (:<, :≀) && shock_first) || (expr.args[1] in (:>, :β‰₯) && !shock_first) + + return shock, is_upper_bound, bound +end + + function mat_mult_kron(A::AbstractArray{T},B::AbstractArray{T},C::AbstractArray{T}; tol::AbstractFloat = eps()) where T <: Real n_rowB = size(B,1) @@ -601,7 +620,7 @@ function match_pattern(strings::Union{Set,Vector}, pattern::Regex) return filter(r -> match(pattern, string(r)) !== nothing, strings) end -function parse_occasionally_binding_constraints(equations_block; max_obc_shift::Int = 20) +function parse_occasionally_binding_constraints(equations_block; max_obc_shift::Int = 10) eqs = [] condition_list = [] @@ -646,7 +665,7 @@ function parse_occasionally_binding_constraints(equations_block; max_obc_shift:: x.args[1] ∈ obc_shocks ? begin obc_shock = intersect([x.args[1]], obc_shocks)[1] - obc_shifts = [Expr(:ref,Meta.parse(string(obc_shock) * "α΅’α΅‡αΆœβ½β»"*super(string(i))*"⁾"),:(x-$i)) for i in 0:max_obc_shift] + obc_shifts = [Expr(:ref,Meta.parse(string(obc_shock) * "α΅’α΅‡αΆœβ½β»"*super(string(i))*"⁾"),i > 0 ? :(x - $i) : :x) for i in 0:max_obc_shift] Expr(:call,:+, x, obc_shifts...) end : x : @@ -2675,7 +2694,7 @@ function solve!(𝓂::β„³; @assert solved "Could not find stable first order solution." - state_update₁ = function(state::Vector{Float64}, shock::Vector{Float64}) sol_mat * [state[𝓂.timings.past_not_future_and_mixed_idx]; shock] end + state_update₁ = function(state::Vector{<: Union{Real,AffExpr}}, shock::Vector{<: Union{Real,AffExpr}}) sol_mat * [state[𝓂.timings.past_not_future_and_mixed_idx]; shock] end 𝓂.solution.perturbation.first_order = perturbation_solution(sol_mat, state_update₁) 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:riccati, :first_order]) @@ -2694,7 +2713,7 @@ function solve!(𝓂::β„³; @assert converged "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1." - state_updateβ‚‚ = function(state::Vector{Float64}, shock::Vector{Float64}) + state_updateβ‚‚ = function(state::Vector{<: Union{Real,AffExpr}}, shock::Vector{<: Union{Real,AffExpr}}) aug_state = [state[𝓂.timings.past_not_future_and_mixed_idx] 1 shock] @@ -2715,7 +2734,7 @@ function solve!(𝓂::β„³; @assert converged "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1." - state_updateβ‚‚ = function(pruned_states::Vector{Vector{Float64}}, shock::Vector{Float64}) + state_updateβ‚‚ = function(pruned_states::Vector{Vector{<: Union{Real,AffExpr}}}, shock::Vector{<: Union{Real,AffExpr}}) aug_state₁ = [pruned_states[1][𝓂.timings.past_not_future_and_mixed_idx]; 1; shock] aug_stateβ‚‚ = [pruned_states[2][𝓂.timings.past_not_future_and_mixed_idx]; 0; zero(shock)] @@ -2736,7 +2755,7 @@ function solve!(𝓂::β„³; @assert converged "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1." - state_update₃ = function(state::Vector{Float64}, shock::Vector{Float64}) + state_update₃ = function(state::Vector{<: Union{Real,AffExpr}}, shock::Vector{<: Union{Real,AffExpr}}) aug_state = [state[𝓂.timings.past_not_future_and_mixed_idx] 1 shock] @@ -2754,7 +2773,7 @@ function solve!(𝓂::β„³; @assert converged "Solution does not have a stochastic steady state. Try reducing shock sizes by multiplying them with a number < 1." - state_update₃ = function(pruned_states::Vector{Vector{Float64}}, shock::Vector{Float64}) + state_update₃ = function(pruned_states::Vector{Vector{<: Union{Real,AffExpr}}}, shock::Vector{<: Union{Real,AffExpr}}) aug_state₁ = [pruned_states[1][𝓂.timings.past_not_future_and_mixed_idx]; 1; shock] aug_state₁̂ = [pruned_states[1][𝓂.timings.past_not_future_and_mixed_idx]; 0; shock] aug_stateβ‚‚ = [pruned_states[2][𝓂.timings.past_not_future_and_mixed_idx]; 0; zero(shock)] @@ -2786,7 +2805,7 @@ function solve!(𝓂::β„³; sol_mat, converged = calculate_quadratic_iteration_solution(βˆ‡β‚; T = 𝓂.timings) - state_updateβ‚β‚œ = function(state::Vector{Float64}, shock::Vector{Float64}) sol_mat * [state[𝓂.timings.past_not_future_and_mixed_idx]; shock] end + state_updateβ‚β‚œ = function(state::Vector{<: Union{Real,AffExpr}}, shock::Vector{<: Union{Real,AffExpr}}) sol_mat * [state[𝓂.timings.past_not_future_and_mixed_idx]; shock] end 𝓂.solution.perturbation.quadratic_iteration = perturbation_solution(sol_mat, state_updateβ‚β‚œ) 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:quadratic_iteration, :binder_pesaran]) @@ -2807,7 +2826,7 @@ function solve!(𝓂::β„³; sol_mat = calculate_linear_time_iteration_solution(βˆ‡β‚; T = 𝓂.timings) - state_updateβ‚β‚œ = function(state::Vector{Float64}, shock::Vector{Float64}) sol_mat * [state[𝓂.timings.past_not_future_and_mixed_idx]; shock] end + state_updateβ‚β‚œ = function(state::Vector{<: Union{Real,AffExpr}}, shock::Vector{<: Union{Real,AffExpr}}) sol_mat * [state[𝓂.timings.past_not_future_and_mixed_idx]; shock] end 𝓂.solution.perturbation.linear_time_iteration = perturbation_solution(sol_mat, state_updateβ‚β‚œ) 𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:linear_time_iteration]) diff --git a/src/macros.jl b/src/macros.jl index ca34d3b1..8dfad5dc 100644 --- a/src/macros.jl +++ b/src/macros.jl @@ -109,6 +109,14 @@ macro model(𝓂,ex...) model_ex, condition_list = parse_occasionally_binding_constraints(model_ex) + obc_shock_bounds = Tuple{Symbol, Bool, Float64}[] + + for c in condition_list + if c isa Expr + push!(obc_shock_bounds, parse_obc_shock_bounds(c)) + end + end + # write down dynamic equations and add auxilliary variables for leads and lags > 1 for (i,arg) in enumerate(model_ex.args) if isa(arg,Expr) @@ -849,6 +857,8 @@ macro model(𝓂,ex...) $T, + $obc_shock_bounds, + solution( perturbation( perturbation_solution(SparseMatrixCSC{Float64, Int64}(β„’.I,0,0), x->x), perturbation_solution(SparseMatrixCSC{Float64, Int64}(β„’.I,0,0), x->x), diff --git a/src/structures.jl b/src/structures.jl index 0e566a35..20dcb8c6 100644 --- a/src/structures.jl +++ b/src/structures.jl @@ -349,6 +349,9 @@ mutable struct β„³ model_third_order_derivatives::Vector{Function} timings::timings + + obc_shock_bounds::Vector{Tuple{Symbol, Bool, Float64}} + solution::solution # symbolics::symbolics diff --git a/test/amend_parser.jl b/test/amend_parser.jl index 1c77f20e..83a1b60d 100644 --- a/test/amend_parser.jl +++ b/test/amend_parser.jl @@ -26,11 +26,13 @@ end get_irf(testmax)(:r,:,:Ο΅αΆ»α΅’α΅‡αΆœβ½β»ΒΉβΎ) get_solution(testmax) -using Optimization, Ipopt, OptimizationMOI, OptimizationOptimJL +# using Optimization, Ipopt, OptimizationMOI, OptimizationOptimJL, LineSearches, OptimizationNLopt, KNITRO import MacroTools: postwalk, unblock import MacroModelling: parse_for_loops, convert_to_ss_equation, simplify,create_symbols_eqs!,remove_redundant_SS_vars!,get_symbols, parse_occasionally_binding_constraints, parse_algorithm_to_state_update import DataStructures: CircularBuffer import Subscripts: super, sub +import LinearAlgebra as β„’ +using BenchmarkTools, JuMP, StatusSwitchingQP 𝓂 = testmax algorithm = :first_order @@ -38,21 +40,267 @@ state_update, pruning = parse_algorithm_to_state_update(algorithm, 𝓂) T = 𝓂.timings periods = 40 -initial_state = zeros(Real, T.nVars) -Y = zeros(Real,T.nVars,periods,1) -T.exo -shock_history = zeros(Real,T.nExo,periods) -shock_history[1,1] = -3 +initial_state = zeros(T.nVars) +# Y = zeros(Real,T.nVars,periods,1) +# T.exo +obc_shocks = [i[1] for i in 𝓂.obc_shock_bounds] + +shock_history = zeros(T.nExo,periods) +shock_history[1,1] = -4 reference_steady_state, solution_error = 𝓂.solution.outdated_NSSS ? 𝓂.SS_solve_func(𝓂.parameter_values, 𝓂, verbose) : (copy(𝓂.solution.non_stochastic_steady_state), eps()) -shock_history[16,1] +# shock_history[16,1] + +obc_shock_idx = contains.(string.(T.exo),"α΅’α΅‡αΆœ") + +function calc_state(x::Vector{S}) where S + Y = zeros(AffExpr,T.nVars,periods,1) + # println(collect(x)) + shock_hst = AffExpr.(copy(shock_history[:,1])) + shock_hst[obc_shock_idx] .= x + + Y[:,1,1] = state_update(initial_state, shock_hst) + + for t in 1:periods-1 + Y[:,t+1,1] = state_update(Y[:,t,1],shock_history[:,t+1]) + end + + Y .+= reference_steady_state[1:T.nVars] + + return Y[4,:,:] +end + +# get bound on shocks +# get minmax condition + +# fill(a[3],sum(obc_shock_idx)Γ·length(𝓂.obc_shock_bounds)) + +model = Model(StatusSwitchingQP.Optimizer) +set_silent(model) + +periods_per_shock = sum(obc_shock_idx)Γ·length(𝓂.obc_shock_bounds) +num_shocks = length(𝓂.obc_shock_bounds) + +# Create the variables over the full set of indices first. +@variable(model, x[1:num_shocks*periods_per_shock]) + +# Now loop through obc_shock_bounds to set the bounds on these variables. +for (idx, v) in enumerate(𝓂.obc_shock_bounds) + is_upper_bound = v[2] + bound = v[3] + idxs = (idx - 1) * periods_per_shock + 1:idx * periods_per_shock + if is_upper_bound + set_upper_bound.(x[idxs], bound) + else + set_lower_bound.(x[idxs], bound) + end +end + +@objective(model, Min, x' * β„’.I * x) +@constraint(model, calc_state(x) .>= 0) +JuMP.optimize!(model) + + +@profview for i in 1:100 JuMP.optimize!(model) end + +value.(x) + +# ≀β‰₯ +:(a-b-c)|>dump +testmax.dyn_equations[3]|>dump + +eq = testmax.dyn_equations[3] +[i for (i,c) in enumerate(condition_list) if c isa Expr] + +import MacroModelling: get_symbols, match_pattern + +#check if equation contains maxmin + +expression_for_placeholder = [] +ismax = [] +eqs = postwalk(x -> + x isa Expr ? + x.head == :call ? + x.args[1] ∈ [:max,:min] ? + begin + if length(intersect(get_symbols(x.args[2]),testmax.var)) == 0 + push!(expression_for_placeholder,x.args[2]) + else + push!(expression_for_placeholder,x.args[3]) + end + if x.args[1] == :max + push!(ismax,true) + :max_placeholder + else + push!(ismax,false) + :min_placeholder + end + end : + x : + x : + x, + eq) + +# solve for the placeholder +using SymPyPythonCall + +@syms max_placeholder, min_placeholder, rβ‚β‚€β‚Ž + +eq_to_solve = eval(eqs) + +if ismax[1] + sol = solve(eq_to_solve, max_placeholder) +else + sol = solve(eq_to_solve, min_placeholder) +end + + +final_expr = Expr(:call, :-, ismax[1] ? :($(expression_for_placeholder[1])) : :(-$(expression_for_placeholder[1])), Meta.parse(string(sol[1]))) + + + + + + +future_varss = collect(reduce(union,match_pattern.(get_symbols.(𝓂.dyn_equations),r"β‚β‚β‚Ž$"))) +present_varss = collect(reduce(union,match_pattern.(get_symbols.(𝓂.dyn_equations),r"β‚β‚€β‚Ž$"))) +past_varss = collect(reduce(union,match_pattern.(get_symbols.(𝓂.dyn_equations),r"β‚β‚‹β‚β‚Ž$"))) +shock_varss = collect(reduce(union,match_pattern.(get_symbols.(𝓂.dyn_equations),r"β‚β‚“β‚Ž$"))) +ss_varss = collect(reduce(union,match_pattern.(get_symbols.(𝓂.dyn_equations),r"β‚β‚›β‚›β‚Ž$"))) + +sort!(future_varss ,by = x->replace(string(x),r"β‚β‚β‚Ž$"=>"")) #sort by name without time index because otherwise eps_zα΄Έβ½β»ΒΉβΎβ‚β‚‹β‚β‚Ž comes before eps_zβ‚β‚‹β‚β‚Ž +sort!(present_varss ,by = x->replace(string(x),r"β‚β‚€β‚Ž$"=>"")) +sort!(past_varss ,by = x->replace(string(x),r"β‚β‚‹β‚β‚Ž$"=>"")) +sort!(shock_varss ,by = x->replace(string(x),r"β‚β‚“β‚Ž$"=>"")) +sort!(ss_varss ,by = x->replace(string(x),r"β‚β‚›β‚›β‚Ž$"=>"")) + +steady_state = [] +for (i, var) in enumerate(ss_varss) + push!(steady_state,:($var = XΜ„[$i])) + # ii += 1 +end + +ii = 1 + +alll = [] +for var in future_varss + push!(alll,:($var = X[$ii])) + ii += 1 +end + +for var in present_varss + push!(alll,:($var = X[$ii])) + ii += 1 +end + +for var in past_varss + push!(alll,:($var = X[$ii])) + ii += 1 +end + +for var in shock_varss + push!(alll,:($var = X[$ii])) + ii += 1 +end + + +# paras = [] +# push!(paras,:((;$(vcat(𝓂.parameters,𝓂.calibration_equations_parameters)...)) = params)) + +paras = [] +for (i, parss) in enumerate(vcat(𝓂.parameters,𝓂.calibration_equations_parameters)) + push!(paras,:($parss = params[$i])) +end + + + + + +mod_func3 = :(function model_jacobian(X::Vector, params::Vector{Real}, XΜ„::Vector) +$(alll...) +$(paras...) +$(𝓂.calibration_equations_no_var...) +$(steady_state...) +$final_expr +end) + + + +intersect(get_symbols(:rΜ„),testmax.var) +# replace maxmin with the side containing parameters only + + +check_minmax() +max(rΜ„,rΜ‚[0]) = r[0] + +r[0] = max(rΜ„,rΜ‚[0]) +r[0] - rΜ„ > 0 +:(begin +r[0] = max(rΜ„,rΜ‚[0]) + v[0] +r[0] - rΜ„ - v[0] > 0 +end) + +using MacroTools +MacroTools.postwalk(:(r[0] = max(rΜ„,rΜ‚[0]) + v[0])) + +r[0] = max(rΜ„,rΜ‚[0]) + v[0] +r[0] - rΜ„ - v[0] > 0 + +:(Ο΅αΆ» > 0) |> dump + + +using MacroTools + +function parse_obc_shock_bounds(expr::Expr) + # Determine the order of the shock and bound in the expression + shock_first = isa(expr.args[2], Symbol) + + # Extract the shock and bound from the expression + shock = shock_first ? expr.args[2] : expr.args[3] + bound_expr = shock_first ? expr.args[3] : expr.args[2] + + # Evaluate the bound expression to get a numerical value + bound = eval(bound_expr) |> Float64 + + # Determine whether the bound is a lower or upper bound based on the comparison operator and order + is_upper_bound = (expr.args[1] in (:<, :≀) && shock_first) || (expr.args[1] in (:>, :β‰₯) && !shock_first) + + return shock, is_upper_bound, bound +end + +# Example usage: +expressions = [:(Ο΅αΆ» > 0), :(1/2 < Ο΅αΆ»), :(Ο΅αΆ» > 1/3+1), :(Ο΅αΆ» ≀ -1), :(1 β‰₯ Ο΅αΆ»)] +parsed_expressions = [parse_expression(expr) for expr in expressions] + + sol_mat = 𝓂.solution.perturbation.first_order.solution_matrix + function bound_violations(x,p) - shock_history[2:end,1] = x + # shock_history[2:end,1] = x + + # Y[:,1,1] = sol_mat * [initial_state[𝓂.timings.past_not_future_and_mixed_idx]; shock_history[:,1]] #state_update(initial_state,shock_history[:,1]) + + # for t in 1:periods-1 + # Y[:,t+1,1] = sol_mat * [Y[𝓂.timings.past_not_future_and_mixed_idx,t,1]; shock_history[:,t+1]] #state_update(Y[:,t,1],shock_history[:,t+1]) + # end + + # Y .+= reference_steady_state[1:T.nVars] + + # target = sum(abs,Y[4,:,:] .- max.(0,Y[5,:,:])) + sum(abs2,x) + return sum(x.^2)#target +end + +# bound_violation(4*zeros(21),[]) - Y[:,1,1] = sol_mat * [initial_state[𝓂.timings.past_not_future_and_mixed_idx]; shock_history[:,1]] #state_update(initial_state,shock_history[:,1]) +# bound_violations([zeros(14)...,2,zeros(6)...],[]) + +function shock_constraints(res, x, p) + # shock_history[2:end,1] = x + + Y = zeros(Real,T.nVars,periods,1) + Y[:,1,1] = sol_mat * [initial_state[𝓂.timings.past_not_future_and_mixed_idx]; vcat(shock_history[1,1],x)] #state_update(initial_state,shock_history[:,1]) for t in 1:periods-1 Y[:,t+1,1] = sol_mat * [Y[𝓂.timings.past_not_future_and_mixed_idx,t,1]; shock_history[:,t+1]] #state_update(Y[:,t,1],shock_history[:,t+1]) @@ -60,23 +308,204 @@ function bound_violations(x,p) Y .+= reference_steady_state[1:T.nVars] - target = sum(abs,Y[4,:,:] .- max.(0,Y[5,:,:])) + sum(abs2,x) - return target + res .= vcat(x, Y[4,:,:]) end -bound_violation(4*zeros(21),[]) +# T.exo -bound_violation([zeros(14)...,2,zeros(6)...],[]) +x0 = zeros(T.nExo-1) +# x0[indexin([:Ο΅αΆ»α΅’α΅‡αΆœβ½β»β°βΎ], T.exo)[1] - 1] = 2 -T.exo[16] -function shock_constraints(res, x, p) - res .= x +# bound_violations(x0*.2,()) +# shock_constraints(zeros(T.nExo-1 + periods), x0, []) + +# optprob = OptimizationFunction(bound_violations, Optimization.AutoForwardDiff(), cons = shock_constraints) +# prob = OptimizationProblem(optprob, x0, (), lcons = zeros(T.nExo-1 + periods), ucons = fill(Inf, T.nExo-1 + periods)) + +# # vcat(zeros(14),3,zeros(6 + periods))[16] +# # T.exo[16] +# # shock_constraints(zeros(21 + periods),vcat(zeros(14),2,zeros(6 + periods)),[]) + +# sol = solve(prob, Ipopt.Optimizer()) + +# sol = solve(prob, MadNLP.Optimizer()) + +# sol = solve(prob, AmplNLWriter.Optimizer("couenne")) + + + + +function calc_state(x::Vector{S}) where S + Y = zeros(AffExpr,T.nVars,periods,1) + + Y[:,1,1] = state_update(initial_state, vcat(shock_history[1,1],x)) + + for t in 1:periods-1 + Y[:,t+1,1] = state_update(Y[:,t,1],shock_history[:,t+1]) + end + + Y .+= reference_steady_state[1:T.nVars] + return Y[4,:,:] +end + + +using BenchmarkTools, JuMP, StatusSwitchingQP +# using BenchmarkTools, JuMP, Ipopt, COSMO, Clarabel, DAQP, HiGHS, MadNLP, OSQP, SCS, StatusSwitchingQP, Hypatia +# import MultiObjectiveAlgorithms as MOA +# model = Model(Ipopt.Optimizer) + +import LinearAlgebra as β„’ + + + +model = Model(StatusSwitchingQP.Optimizer) +set_silent(model) +@variable(model, x[i = 1:T.nExo - 1] >= 0) +@objective(model, Min, x' * β„’.I * x) +@constraint(model, calc_state(x) .>= 0) +JuMP.optimize!(model) + +value.(x) + +@benchmark begin + model = Model(COSMO.Optimizer) + set_silent(model) + @variable(model, x[1:11] >= 0) + @objective(model, Min, x'*Q*x) + @constraint(model, calc_state(x) .>= 0) + JuMP.optimize!(model) end -x0 = zeros(21) .+ .1 -optprob = OptimizationFunction(bound_violations, Optimization.AutoForwardDiff(), cons = shock_constraints) -prob = OptimizationProblem(optprob, x0, (), lcons = zeros(21), ucons = fill(Inf,21)) +@benchmark begin + model = Model(MadNLP.Optimizer) + set_silent(model) + @variable(model, x[1:11] >= 0) + @objective(model, Min, x'*Q*x) + @constraint(model, calc_state(x) .>= 0) + JuMP.optimize!(model) +end + +# favorite +@benchmark begin + model = Model(StatusSwitchingQP.Optimizer) + set_silent(model) + @variable(model, x[1:11] >= 0) + @objective(model, Min, x'*Q*x) + @constraint(model, calc_state(x) .>= 0) + JuMP.optimize!(model) +end + +value.(x) + +@benchmark begin + model = Model(Clarabel.Optimizer) + set_silent(model) + @variable(model, x[1:11] >= 0) + @objective(model, Min, x'*Q*x) + @constraint(model, calc_state(x) .>= 0) + JuMP.optimize!(model) +end + + +@benchmark begin + model = Model(Hypatia.Optimizer) + set_silent(model) + @variable(model, x[1:11] >= 0) + @objective(model, Min, x'*Q*x) + @constraint(model, calc_state(x) .>= 0) + JuMP.optimize!(model) +end + + + + + +@benchmark begin + model = Model(SCS.Optimizer) + set_silent(model) + @variable(model, x[1:11] >= 0) + @objective(model, Min, x'*Q*x) + @constraint(model, calc_state(x) .>= 0) + JuMP.optimize!(model) +end + + +@benchmark begin + model = Model(OSQP.Optimizer) + set_silent(model) + @variable(model, x[1:11] >= 0) + @objective(model, Min, x'*Q*x) + @constraint(model, calc_state(x) .>= 0) + JuMP.optimize!(model) +end + + +@benchmark begin + model = Model(DAQP.Optimizer) + set_silent(model) + @variable(model, x[1:11] >= 0) + @objective(model, Min, x'*Q*x) + @constraint(model, calc_state(x) .>= 0) + JuMP.optimize!(model) +end + + +@benchmark begin + model = Model(HiGHS.Optimizer) + set_silent(model) + @variable(model, x[1:11] >= 0) + @objective(model, Min, x'*Q*x) + @constraint(model, calc_state(x) .>= 0) + JuMP.optimize!(model) +end + + + +@benchmark begin + model = Model(Ipopt.Optimizer) + set_silent(model) + @variable(model, x[1:11] >= 0) + @objective(model, Min, x'*Q*x) + @constraint(model, calc_state(x) .>= 0) + JuMP.optimize!(model) +end + + +@benchmark begin + model = Model(ProxSDP.Optimizer) + set_silent(model) + @variable(model, x[1:11] >= 0) + @objective(model, Min, x'*Q*x) + @constraint(model, calc_state(x) .>= 0) + JuMP.optimize!(model) +end + + + + +solution_summary(model) + +value.(x) + + + +# sol = solve(pro, SCIP.Optimizer()) +# sol = solve(prob, Pavito.Optimizer()) +# sol = solve(prob, Juniper.Optimizer()) +# sol = solve(prob, COSMO.Optimizer()) +# sol = solve(prob, BARON.Optimizer()) +# sol = solve(prob, knitro.Optimizer()) +# sol = solve(prob, highs.Optimizer()) +# sol = solve(prob, EAGO.Optimizer()) + +bound_violations(sol*.9, ()) +shock_constraints(zeros(T.nExo-1 + periods), sol, []) + +sol = solve(prob, NelderMead()) + +sol = solve(prob, LBFGS(linesearch = LineSearches.BackTracking(order=3))) + sol = solve(prob, IPNewton()) @@ -162,6 +591,15 @@ eq = :(begin end) model_ex, condition_list = parse_occasionally_binding_constraints(eq) + +obc_shock_bounds = Tuple{Symbol, Bool, Float64}[] + +for c in condition_list + if c isa Expr + push!(obc_shock_bounds, parse_obc_shock_bounds(c)) + end +end + obc_shocks = Symbol[] for a in condition_list