# Condition for particular index in array to be maximum

Hi,

I am new to Julia and trying to solve a fitting problem. I have a transition matrix that I am trying to fit. I have got some results but they are non-intutive as the probability of being in state i,i should be the greatest. How can I put such a constraint in JuMP?

The working code (not tested) is as follows:

``````function calibrate_trans_rate(filepath_pd)

# Get PD values
pd_table = zeros(Float64, 30, 7)
open(filepath_pd, "r") do fp
for row in 1:30
pd_table[row, :] .= parse.(Float64, split(line, ","))
end
end

# Set up solver
model = Model(Ipopt.Optimizer)  # need to decide what solver to use

# Set variable to be transition values
@variable(model, 0 <= x[1:8, 1:8] <= 1)
@variable(model, y[1:8, 1:8, 1:31])	# This is needed as JuMP only supports matrix operators up to quadratic terms. It doesn’t support cubic.

# Constraints on transition rates
@constraints(model, begin
default_col[i = 1:7], x[i, 8] == pd_table[1, i]
default_row[i = 1:7], x[8, i] == 0
default_state, x[8, 8] == 1
sum_of_rows[i = 1:8], sum(x[i, :]) == 1
end)

# Constraints on y {this is just a hack to get cubic terms for JuMP}
@constraints(model, begin
first_y_cons, y[:, :, 1] .== x  # Or maybe set equal to x
[t = 1:30], x * y[:, :, t] .== y[:, :, t + 1]
end)

# Define statistic to minimize
@NLexpression(model, tot_sum_of_sq_diff, sum((y[k,8,time] - pd_table[1, k])^2 for k=1:7, time=1:30))

@NLobjective(model, Min, tot_sum_of_sq_diff)

JuMP.optimize!(model)

obj_value = JuMP.objective_value(model)
trans_matrix = zeros(8, 8)
for row=1:8, col=1:8
trans_matrix[row, col] = JuMP.value(x[row,col])
end

return trans_matrix
end
``````

I tried adding a line in @constraints as follows:

``````curr_state[i = 1:8], x[i,i] == maximum(x[i, :])
``````

This gives me the following error:

``````MethodError: no method matching isless(::JuMP.VariableRef, ::JuMP.VariableRef)

Closest candidates are:

isless(!Matched::Missing, ::Any) at missing.jl:87

isless(::Any, !Matched::Missing) at missing.jl:88

max(::JuMP.VariableRef, ::JuMP.VariableRef)@operators.jl:417
_mapreduce@reduce.jl:408[inlined]
_mapreduce_dim@reducedim.jl:318[inlined]
#mapreduce#620@reducedim.jl:310[inlined]
mapreduce@reducedim.jl:310[inlined]
_maximum@reducedim.jl:727[inlined]
_maximum@reducedim.jl:726[inlined]
#maximum#631@reducedim.jl:722[inlined]
maximum@reducedim.jl:722[inlined]
macro expansion@rewrite.jl:227[inlined]
macro expansion@macros.jl:440[inlined]
(::Main.workspace192.var"#7#15"{JuMP.Model,Array{JuMP.VariableRef,2}})(::Int64)@macro.jl:183
#26@container.jl:70[inlined]
iterate@generator.jl:47[inlined]
collect(::Base.Generator{JuMP.Containers.VectorizedProductIterator{Tuple{Base.OneTo{Int64}}},JuMP.Containers.var"#26#27"{Main.workspace192.var"#7#15"{JuMP.Model,Array{JuMP.VariableRef,2}}}})@array.jl:686
map@abstractarray.jl:2188[inlined]
container@container.jl:70[inlined]
container@container.jl:65[inlined]
macro expansion@macros.jl:91[inlined]
calibrate_trans_rate(::String)@Other: 25
top-level scope@Local: 1
``````

Any ideas?

The problem that I am trying to solve is here:

Ignore the question - the error was in this line and pd_table should be:

``````@NLexpression(model, tot_sum_of_sq_diff, sum((y[k,8,time] - pd_table[time, k])^2 for k=1:7, time=1:30))
``````
1 Like

You probably want

``````@constraint(model, [i = 1:8, j = 1:8], x[i, i] >= x[i, j])
``````