Hello,
I am trying to solve an MCP that requires performing operations across some indices.
Here is a MWE:
using JuMP
using Random
using LinearAlgebra
using BenchmarkTools
using JuMP
import PATHSolver
PATHSolver.c_api_License_SetString("2830898829&Courtesy&&&USR&45321&5_1_2021&1000&PATH&GEN&31_12_2025&0_0_0&6000&0_0")
struct MyArrays
denom_i::Vector{Float64}
v::Vector{Float64}
u::Vector{Float64}
end
function make_arrays(S, J)
MyArrays(zeros(S), zeros(S), zeros(J))
end
function example_one_go(J::Int, S::Int, A::Matrix, Guess_y::Vector{Float64}, z::Vector{Float64}, C::Matrix, Guess::Matrix, myarrays::MyArrays)
model = Model(PATHSolver.Optimizer)
set_silent(model)
@variable(model, x[i = 1:J, j = 1:J] >= 0, start = Guess[i, j])
y_fixed_1 = 1.0
@variable(model, y[i = 2:J] >= 0, start = Guess_y[i])
last_x = Vector{Float64}()
function f(x_i, x_j, x...)
i = round(Int, x_i)
j = round(Int, x_j)
new_x_vec = collect(x)
if last_x === nothing || new_x_vec != last_x
X = reshape(new_x_vec, J, J)
denom_i = myarrays.denom_i
@views mul!(denom_i, A, X[i, :])
last_x = x
end
return sum(A[s, j] / denom_i[s] for s in 1:S) * (1 / S)
end
function ∇f(g, x_i, x_j, x...)
local i = round(Int, x_i)
local j = round(Int, x_j)
new_x_vec = collect(x)
if last_x === nothing || new_x_vec != last_x
X = reshape(new_x_vec, J, J)
fill!(g, 0.0)
denom_i = myarrays.denom_i
@views mul!(denom_i, A, X[i, :])
v = myarrays.v
u = myarrays.u
@inbounds @simd for s in 1:S
v[s] = A[s, j] / (denom_i[s]^2)
end
@views mul!(u, transpose(A), v)
for q in 1:J
local col = 2 + (q - 1) * J + i
g[col] = -(1 / S) * u[q]
end
last_x = x
end
return
end
get_y(i) = (i == 1) ? y_fixed_1 : y[i]
@operator(model, op_f, 2 + J * J, f, ∇f)
@constraint(model, comp1[i in 1:J, j in 1:J], -get_y(i) * op_f(i, j, x...) + get_y(j) * C[i, j] / z[j] ⟂ x[i, j])
@constraint(model, comp2[i in 2:J], -sum(x[j, i] * C[i, j] for j in 1:J) + z[i] ⟂ y[i])
JuMP.optimize!(model)
Y = vcat(1.0, collect(value.(y)))
return value.(x), Y, ((Y ./ z) .* C) .* value.(x)
end
# Example inputs to run the function
J = 10
S = 100
Guess_y = ones(J)
z = normalize(rand(J), 1) .+ 0.1 # Normalize z to have unit sum
C = ones(J, J)
for i in 1:J
for j in 1:J
if i != j
C[i, j] = 1.25
end
end
end
Guess = ones(J, J)
A = rand(S,J) .+ 0.1 # Normalize and shift values to avoid zeros
myarrays = make_arrays(S, J)
# Running the function
result_x, result_Y, transformed_matrix = example_one_go(J, S, A, Guess_y, z, C, Guess, myarrays)
I have a couple of questions:
-
Are there any performance gains to be achieved in this code? In my true application, both J and S are in the thousands.
-
Is this code amenable to using other MCP solvers such as GAMS or KNitro? Or the syntax/JuMP structure does not translate well?
-
Is there any way that is useful for performance in which instead of looping over i,j, I just do everything once filling up a large (possible SParse) matrix, and then call that matrix to evaluate the constraints?