Condition for particular index in array to be maximum

Specific Domains Optimization (Mathematical)
,
1

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
			line = readline(fp)
			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:

2

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
3

Related to Array variable - matrix multplication in JuMP - #3 by shashankchahar

You probably want

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