Help with a Constraint

Specific Domains Optimization (Mathematical)
1

Hello,

I need some help to write a constraint for my problem.
This constraint has the following matrix:

using JuMP, CouenneNL

Nt = 2;
L = 0;
C = 0;
h = 0;
hmax = 5;
Y = zeros(Nt,Nt,hmax);
b = 2;
v = 10;

for h = 3:2:hmax
  for L = 1:Nt
    for C = 1:Nt
      Y[L,C,h] = b*h;
      b = b+1;
    end
  end
  b = 2;
end

But in the constraint, this matrix has to be modified a little bit. I need to insert binary variables (say x(i=1:Nt)) in the main diagonal (this binary variable means that if its equal 1 i’ll install certain equipment, and if it’s 0 i won’t).

To be more clear, it would be something like this:

Y = [2*h+xi*v     3*h
       4h      5*h+xi*v]

But i can only use this variable in the optimizaton problem, i can’t use it before.

So, we got h = 3:2:hmax and my constraint is something like:

m = Model()
@constraint(m,Y[:,:,h]*U[:,h] - I[:,h] == 0)

Where:

U is an (Nt X hmax) array, the solution of my problem
I is  an (Nt X hmax) data array that was previously determined

How can i insert the binary variables in that constraint?

I don’t know if i made myself very clear. If you have any questions, i’ll gladly answer them.

Thanks in advance!

2

You can add a vector of Nt binary variables as:

Nt = 2
model = Model()
@variable(model, x[i=1:Nt], Bin)

I can’t help with the rest of your problem since I don’t know what xi, U, and I are.

3

Hello odow, thanks for your reply.

I’m sorry, that was a typo, it’s not xi.

It’s x[i], the vector of binary variables, like you just defined. The question is, how can i add those variables in the main diagonal of my Y matrix in this constraint?

@constraint(m,Y[:,:,h]*U[:,h] - I[:,h] == 0)

U is the vector of the variables i wish to determinate
I is a vector of data input i already have.

And i have another question… supose i have this variable hs:

3 <= hs <= 25

But hs can only assume odd values between those limites, i.e, [3, 5, 7, ..., 25]. How can this be defined?

4

It looks like you might be new to JuMP and formulating problems as constrained optimization problems. You may want to have a read of the JuMP docs: Quick Start Guide — JuMP -- Julia for Mathematical Optimization 0.18 documentation

There are also some example notebooks: https://github.com/JuliaOpt/juliaopt-notebooks/tree/3110eddaf5effdecfee687739295bea05731ba33/notebooks

But hs can only assume odd values between those limites, i.e, [3, 5, 7, ..., 25] . How can this be defined?

You could add a dummy variable and constraint like this:

model = Model()
@variable(model, 3 <= hs <= 25, Int)
@variable(model, hs_dummy, Int)
@constraint(model, hs == 2 * hs_dummy + 1)
1 Like
5

Allright, thanks a lot for your help! Will read those links now!

I was trying work in some variables, but i got some errors:

hsmin = 3*ones(Nt,1);
hsmax = 25*ones(Nt,1);

m = Model()

@variable(m, hsmin[i] <= hs[i=1:Nt] <= hsmax[i], Int)

I got this error:

ERROR: LoadError: The operators <=, >=, and == can only be used to specify scalar constraints. If you are trying to add a vectorized constraint, use the element-wise
dot comparison operators (.<=, .>=, or .==) instead

Then i tried the element-wise comparison, like this:

hsmin = 3*ones(Nt,1);
hsmax = 25*ones(Nt,1);

m = Model()

@variable(m, hsmin[i] .<= hs[i=1:Nt] .<= hsmax[i], Int)

And i got this error:

ERROR: LoadError: In @variable(m,hsmin[i] .<= hs[i = 1:Nt] .<= hsmax[i],Int): Use the form lb <= ... <= ub.

6

This works for me:

Nt = 3
hsmin = 3*ones(Nt,1)
hsmax = 25*ones(Nt,1)

m = Model()

@variable(m, hsmin[i] <= hs[i=1:Nt] <= hsmax[i], Int)
7

Thats weird… i tried today, and it worked. Well, moving on…

I’ve added some other variables/constraints to my problem:

using JuMP

Nt = 2;
L = 0;
C = 0;
h = 0;
hmax = 5;
Yh = zeros(Nt,Nt,hmax);
b = 2;

for h = 3:2:hmax
  for L = 1:Nt
    for C = 1:Nt
      Yh[L,C,h] = b*h;
      b = b+1;
    end
  end
  b = 2;
end

Ih = [2 1 3 4 5; 7 5 6.5 3 10];
hsmin = 3*ones(Nt,1);
hsmax = 25*ones(Nt,1);

m = Model(solver=CouenneNLSolver())

@variable(m, hsmin[i] <= hs[i=1:Nt] <= hsmax[i], Int)
@variable(m, hs_dummy[i=1:Nt], Int)
@variable(m, Vh[i = 1:Nt,h = 3:2:hmax])

@NLobjective(m, Min, sum(10*hs[i] for i = 1:Nt))
@constraint(m, hs[i=1:Nt] .== 2 * hs_dummy[i] + 1)
@constraint(m, con[i = 1:Nt,h = 3:2:hmax], Yh[i,i,h]*Vh[i,h] == Ih[i,h]) # New Constraint Added
#print(m)

status = solve(m)

println("Objective value: ", getobjectivevalue(m))
println("hs = ", getvalue(hs))
println("Vh = ", getvalue(Vh))

I get as result, this:

Objective value: 60.0
hs = [3.0, 3.0]
Vh = Vh: 2 dimensions:
[1,:]
  [1,3] = 0.5
  [1,5] = 0.5
[2,:]
  [2,3] = 0.43333333333333335
  [2,5] = 0.4

But, if i type Vh[1,3] in Julia, i get the following:

julia> Vh[1,3]
Vh[1,3]

How can i get the actual value of the Vh[1,3], i.e, 0.5? I’m asking this, because later i’ll have to solve a similar but much larger problem, and i don’t want to print all those values of Vh when i solve it, but would be nice to know any of them if i wish so.

8

You have two options:

getvalue(Vh[1, 3])

or

Vh_solution = getvalue(Vh)
Vh_solution[1, 3]
1 Like
9

Hello again, supose i have these parameters and variables:

L = 0;
C = 0;
i = 0;
Nt = 2;
hmax = 5;
@variable(model, a[i=1:Nt], Bin)
Yh = zeros(Complex128,Nt,Nt,hmax);
y1 = [0+0im    0+0im    0.05 - 0.15000im    0+0im    0.01923-0.09615im; 
      0+0im    0+0im    0.07692-0.11538im   0+0im    0.03448-0.0862im];
y2 = [0+0im    0+0im    0+6im    0+0im    0+10im; 
      0+0im    0+0im    0+6im    0+0im    0+10im];

Is there anyway i can write this loop as a constraint?

for h = 3:2:hmax
Yaux = zeros(Complex128,Nt,Nt,h);
  for L = 1:Nt
    for C = 1:Nt
      for i = 1:Nt
        if CONDITION 1
            if CONDITION 2
                 Yaux[L,C,h] = Yaux[L,C,h]+y1[i,h]+a[i]*y2[i,h];
            end
            if CONDITION 3
                Yaux[L,C,h] = Yaux[L,C,h]+y1[i,h];
            end
        end
        if CONDITION 4
                 Yaux[L,C,h] = Yaux[L,C,h]-y1[i,h];
                 Yaux[C,L,h] = Yaux[L,C,h];
        end
    Yh[L,C,h] = Yaux[L,C,h];
      end
    end
  end
end

I’ve read about the sum() syntax in the quick start guide, but couldn’t make it work properly. It’s not possible? or i’m doing it wrong?

10

If you want to construct a constraint that way you should consider using the macro @expression (https://jump.readthedocs.io/en/latest/refexpr.html) to build a affine expression and then you can use it to create a linear constraint.
I would also recommend using the simplest example you can here :slight_smile:

11

Ok, i will look into it, thanks for the reply!

I tried adding another constraint to my problem. Like this:

using JuMP, CouenneNL

Nt = 2;
L = 0;
C = 0;
h = 0;
hmax = 5;
Yh = zeros(Nt,Nt,hmax);
b = 2;

for h = 3:2:hmax
  for L = 1:Nt
    for C = 1:Nt
      Yh[L,C,h] = b*h;
      b = b+1;
    end
  end
  b = 2;
end

Ih = [2 1 3 4 5; 7 5 6.5 3 10];
u = [10;14];
kj = 1;
Qmod = 50;
hsmin = 3*ones(Nt,1);
hsmax = 25*ones(Nt,1);

m = Model(solver=CouenneNLSolver())

@variable(m, hsmin[i] <= hs[i=1:Nt] <= hsmax[i], Int)
@variable(m, hs_dummy[i=1:Nt], Int)
@variable(m, Vh[i = 1:Nt,h = 3:2:hmax])
@variable(m, DHTv[i = 1:Nt])

@NLobjective(m, Min, sum(10*hs[i] for i = 1:Nt))

@constraint(m, hs[i = 1:Nt] .== 2 * hs_dummy[i] + 1)
@constraint(m, con[i = 1:Nt,h = 3:2:hmax], Yh[i,i,h]*Vh[i,h] == Ih[i,h])
@constraint(m, conn[i = 1:Nt, h = 3:2:hmax], DHTv[i] == sum(Vh[i,h]^2)/u[i]^2) # New Constraint added

status = solve(m)

println("Objective value: ", getobjectivevalue(m))
println("hs = ", getvalue(hs))
println("Vh = ", getvalue(Vh))
println("DHTv = ", getvalue(DHTv))

But i’m getting the following error:

ERROR: LoadError: AssertionError: line[1:7] == "Options"
Stacktrace:
 [1] read_sol(::IOStream, ::AmplNLWriter.AmplNLMathProgModel) at C:\Users\Muril\.julia\v0.6\AmplNLWriter\src\AmplNLWriter.jl:666
 [2] read_results(::IOStream, ::AmplNLWriter.AmplNLMathProgModel) at C:\Users\Muril\.julia\v0.6\AmplNLWriter\src\AmplNLWriter.jl:575
 [3] open(::AmplNLWriter.##117#118{AmplNLWriter.AmplNLMathProgModel}, ::String, ::String) at .\iostream.jl:152
 [4] read_results(::AmplNLWriter.AmplNLMathProgModel) at C:\Users\Muril\.julia\v0.6\AmplNLWriter\src\AmplNLWriter.jl:569
 [5] optimize!(::AmplNLWriter.AmplNLMathProgModel) at C:\Users\Muril\.julia\v0.6\AmplNLWriter\src\AmplNLWriter.jl:419
 [6] optimize!(::CouenneNL.CouenneNLNonlinearModel) at C:\Users\Muril\.julia\v0.6\Lazy\src\macros.jl:268
 [7] #solvenlp#165(::Bool, ::Function, ::JuMP.Model, ::JuMP.ProblemTraits) at C:\Users\Muril\.julia\v0.6\JuMP\src\nlp.jl:1271
 [8] (::JuMP.#kw##solvenlp)(::Array{Any,1}, ::JuMP.#solvenlp, ::JuMP.Model, ::JuMP.ProblemTraits) at .\<missing>:0
 [9] #solve#116(::Bool, ::Bool, ::Bool, ::Array{Any,1}, ::Function, ::JuMP.Model) at C:\Users\Muril\.julia\v0.6\JuMP\src\solvers.jl:172
 [10] solve(::JuMP.Model) at C:\Users\Muril\.julia\v0.6\JuMP\src\solvers.jl:150
 [11] include_string(::String, ::String) at .\loading.jl:522
 [12] include_string(::Module, ::String, ::String) at C:\Users\Muril\.julia\v0.6\Compat\src\Compat.jl:88
 [13] (::Atom.##112#116{String,String})() at C:\Users\Muril\.julia\v0.6\Atom\src\eval.jl:109
 [14] withpath(::Atom.##112#116{String,String}, ::String) at C:\Users\Muril\.julia\v0.6\CodeTools\src\utils.jl:30
 [15] withpath(::Function, ::String) at C:\Users\Muril\.julia\v0.6\Atom\src\eval.jl:38
 [16] hideprompt(::Atom.##111#115{String,String}) at C:\Users\Muril\.julia\v0.6\Atom\src\repl.jl:67
 [17] macro expansion at C:\Users\Muril\.julia\v0.6\Atom\src\eval.jl:106 [inlined]
 [18] (::Atom.##110#114{Dict{String,Any}})() at .\task.jl:80
while loading C:\Users\Muril\Desktop\Murilo\UFPE\Dissertação\Julia\Arquivos Feitos\teste_otimização_matriz_3_dimensões.jl, in expression starting on line 56

Does this mean that my problem is now infeasible?

I’m not worried at the moment if the problem is feasible or not, i’m just trying to learn how to write those @variable, @constraint and @objective macros.

Is this written correct? I mean, it would work, if the supposed infeasibility of the problem wasn’t causing the error here?

@constraint(m, conn[i = 1:Nt, h = 3:2:hmax], DHTv[i] == sum(Vh[i,h]^2)/u[i]^2) # New Constraint added

12

Well, the part above is fixed, but i ran into another problem. Heres the updated code:

using JuMP, CouenneNL

Nt = 2;
L = 0;
C = 0;
hmax = 5;
Yh = zeros(Nt,Nt,hmax);
b = 2;
Ih = [2 0 3 0 5 3 9; 7 0 6.5 0 10 4 6];
u = [10;14];
hsmin = 3*ones(Nt,1);
hsmax = 25*ones(Nt,1);

for h = 3:2:hmax
  for L = 1:Nt
    for C = 1:Nt
      Yh[L,C,h] = b*h+2.3;
      b = b+1;
    end
  end
  b = 2;
end

m = Model(solver=CouenneNLSolver())

@variable(m, a[i=1:Nt], Bin) # Not using this variable yet
@variable(m, hsmin[i] <= hs[i=1:Nt] <= hsmax[i], Int)
@variable(m, hs_dummy[i=1:Nt], Int)
@variable(m, Vh[i = 1:Nt,h = 3:2:hmax])
@variable(m, DHTv[i = 1:Nt])
@variable(m, DHIv[i = 1:Nt, h = 3:2:hmax])

@NLobjective(m, Min, sum(10*hs[i] for i = 1:Nt))

@constraint(m, hs .== 2 * hs_dummy + 1)
@constraint(m, con1[i = 1:Nt, h = 3:2:hmax], Yh[i,i,h]*Vh[i,h] == Ih[i,h])
@NLconstraint(m, con2[i = 1:Nt], DHTv[i] == sqrt(sum(Vh[i,h]^2 for h = 3:2:hmax))/u[i]^2)
@constraint(m, con3[i = 1:Nt, h = 3:2:hmax], DHIv[i,h] == Vh[i,h]/u[i]) # This constraint is causing problem

status = solve(m)

println("Objective value: ", getobjectivevalue(m))
println("hs = ", getvalue(hs))
println("Vh = ", getvalue(Vh))
println("DHTv = ", getvalue(100*DHTv))
println("DHIv = ", getvalue(DHIv))

If we remove the last constraint here, the code works fine. But after adding it, i got the following error:

ANALYSIS TEST: Error evaluating constraint 1: can't evaluate sqrt'(0).
WARNING: Not solved to optimality, status: Error
Objective value: NaN
WARNING: Variable value not defined for component of hs. Check that the model was properly solved.
hs = [NaN, NaN]
WARNING: Variable value not defined for Vh. Check that the model was properly solved.
Vh = Vh: 2 dimensions:
[1,:]
  [1,3] = NaN
  [1,5] = NaN
[2,:]
  [2,3] = NaN
  [2,5] = NaN
WARNING: Variable value not defined for DHTv[1]. Check that the model was properly solved.
WARNING: Variable value not defined for DHTv[2]. Check that the model was properly solved.
DHTv = [NaN, NaN]
WARNING: Variable value not defined for DHIv. Check that the model was properly solved.
DHIv = DHIv: 2 dimensions:
[1,:]
  [1,3] = NaN
  [1,5] = NaN
[2,:]
  [2,3] = NaN
  [2,5] = NaN

It says it cant evaluate sqrt(0) in the constraint 1. I don’t know how Couenne check those restrictions, but i assume its not by the order they are defined here? Since the first one doesn’t even have the sqrt(). And it’s a bit odd, it can’t evaluate sqrt(0).

Any ideas how to fix this?

EDIT: I found this in the internet:

By default, the initial values of all the variables are zero. But to get
started, MINOS must compute the gradients of the constraint functions, and
that requires it to compute the derivative of the sqrt function at zero –
sqrt’(0) in the error message – which fails because that derivative of sqrt
at zero is infinite.

This was related to MINOS solver, but i guess the same thing happens in Couenne. Just add a really small value (1e-20 for example) to the sqrt() argument, and it should work.

13

Hello, i got stuck again trying to run my code, it says my problem is infeasible. So, i got two questions, if anyone could help me out.

1st: Imagine i got something like this:

Rmax = 3; # Rows
Cmax = 3; # Columns
hmax = 5; # Third index

And i got these previously calculated matrixs, say A[R = 1:Rmax, C = 1:Cmax, h = 3:2:hmax] and another one, say m[R = 1:Rmax, h = 3:2:hmax].

Now, in my optimization problem, i need to insert some variables in the main diagonal of that matrix, to create some restrictions, so i tried this:

@variable(m, a[R=1:Rmax], Bin)
@variable(m, new_A[R = 1:Rmax, C = 1:Cmax, h = 3:2:hmax])

for h = 3:2:hmax, R = 1:Rmax, C = 1:Cmax
   if R == C
       @NLconstraint(m, new_A[R, C, h] == A[R,C,h]+a[R]*m[R,h]);
   else
       @constraint(m, new_A[R, C, h] == A[R,C,h]);
   end
end

2nd: Now, supose from those new_A matrixes, i need new constraints, just like in LP form A*x = b, so i tried this:

for h = 3:2:hmax
    @constraint(m, new_A[:,:,h]*x[:,h] .== b[:,h])
end

Any ideas what could be the problem here?

14

It looks like you have a model called m and a matrix also called m.

It’s hard to say more without a small reproducible example with the actual error.

15

Actually, the m was just a random letter i put in there trying to explain, i could rename it for anything else.

I don’t know how could i put a small example here, my code is not so small and im using some .txt data and a function i created to do some calculations, so i would need to attach those also.

I’m more concerned if i’m writting those constraints in the right way (syntactically speaking), since i don’t know much about Julia and programming in general.

16

Hello,

I’ve tried everything i could and i still cant make it to work. If anyone has a spare time to help me out, would be awesome.

I can send you (via e-mail or something else) my model i wrote in a .pdf file, so you can understand it better and the .jl file wich i was working on, so you can see whats wrong there.

Just let me know, thanks in advance!