Multiplying random variables will result in an error. This is related to nonlinear programming.

[wrx1990]

My model involves the multiplication of results with random variables and other variables, followed by an error occurring.

MathOptInterface.UnsupportedConstraint{MathOptInterface.ScalarQuadraticFunction{Float64}, MathOptInterface.EqualTo{Float64}}: `MathOptInterface.ScalarQuadraticFunction{Float64}`-in-`MathOptInterface.EqualTo{Float64}` constraint is not supported by the model.

the code like this

using SDDP, HiGHS, Test
function demo()
    model = SDDP.LinearPolicyGraph(;
        stages = 5,
        upper_bound = 30.0,
        sense = :Max,
        optimizer = HiGHS.Optimizer,
    ) do sp, t
        p1=2.0
        @variable(sp, 0 <= x1 , SDDP.State, initial_value = 0.0)
        @variable(sp, 0 <= x2 , SDDP.State, initial_value = 0.0)
        @variable(sp, 0 <= x3 , SDDP.State, initial_value = 10.0)
        @variables(sp, begin
            fx1     
            fx2     
            cash    
        end)
        @variable(sp, p2)
        Ω = [-1.75,0.08,1.93]
        P = [0.25,0.5,0.25]
        SDDP.parameterize(sp, Ω, P) do ω
            return JuMP.fix(p2,  ω+p1*0.8)
        end
        @constraints(sp, begin
            x1.out == x1.in -fx1
            x2.out == x2.in -fx2
            cash == p1*fx1 +p2*fx2
            x3.out == x3.in +cash 
        end)
        if t <5
            @stageobjective(sp,  cash)
        else
            @stageobjective(sp,  cash+x1.out*p1+x2.out*p2)
        end
    end
   SDDP.train(model)
    return
end
demo()

I have reviewed solutions to similar issues in the problem thread, but they do not apply to my specific case. What should I do to resolve this problem?

[odow]

JuMP implements `fixed variables by adding a new decision variable with fixed bounds.

This means that p2 * x is a quadratic function.

HiGHS does not support quadratic constraints.

Instead, you should use set_normalized_coefficient and set_objective_coefficient

I didn't run your code, but something like this should work:

using SDDP, HiGHS, Test
function demo()
    model = SDDP.LinearPolicyGraph(;
        stages = 5,
        upper_bound = 30.0,
        sense = :Max,
        optimizer = HiGHS.Optimizer,
    ) do sp, t
        p1 = 2.0
        p2 = 1.0  # PLACEHOLDER
        @variable(sp, 0 <= x1, SDDP.State, initial_value = 0.0)
        @variable(sp, 0 <= x2, SDDP.State, initial_value = 0.0)
        @variable(sp, 0 <= x3, SDDP.State, initial_value = 10.0)
        @variables(sp, begin
            fx1     
            fx2     
            cash    
        end)
        @constraints(sp, begin
            x1.out == x1.in - fx1
            x2.out == x2.in - fx2
            c_cash, cash == p1 * fx1 + p2 * fx2
            x3.out == x3.in + cash 
        end)
        if t < 5
            @stageobjective(sp,  cash)
        else
            @stageobjective(sp, cash + p1 * x1.out + p2 * x2.out)
        end
        Ω = [-1.75, 0.08, 1.93]
        P = [0.25, 0.5, 0.25]
        SDDP.parameterize(sp, Ω, P) do ω
            p2 = ω + p1 * 0.8
            # Note the `-p2` becuase `fx2` appears on the right-hand side
            set_normalized_coefficient(c_cash, fx2, -p2)
            if t == 5
                set_objective_coefficient(sp, x2.out, p2)
            end
            return
        end
    end
    SDDP.train(model)
    return
end

demo()

[odow]

Any other questions? If not, I will close this issue.

[wrx1990]

Any other questions? If not, I will close this issue.

JuMP implements `fixed variables by adding a new decision variable with fixed bounds.

This means that p2 * x is a quadratic function.

HiGHS does not support quadratic constraints.

Instead, you should use set_normalized_coefficient and set_objective_coefficient

I didn't run your code, but something like this should work:

using SDDP, HiGHS, Test
function demo()
    model = SDDP.LinearPolicyGraph(;
        stages = 5,
        upper_bound = 30.0,
        sense = :Max,
        optimizer = HiGHS.Optimizer,
    ) do sp, t
        p1 = 2.0
        p2 = 1.0  # PLACEHOLDER
        @variable(sp, 0 <= x1, SDDP.State, initial_value = 0.0)
        @variable(sp, 0 <= x2, SDDP.State, initial_value = 0.0)
        @variable(sp, 0 <= x3, SDDP.State, initial_value = 10.0)
        @variables(sp, begin
            fx1     
            fx2     
            cash    
        end)
        @constraints(sp, begin
            x1.out == x1.in - fx1
            x2.out == x2.in - fx2
            c_cash, cash == p1 * fx1 + p2 * fx2
            x3.out == x3.in + cash 
        end)
        if t < 5
            @stageobjective(sp,  cash)
        else
            @stageobjective(sp, cash + p1 * x1.out + p2 * x2.out)
        end
        Ω = [-1.75, 0.08, 1.93]
        P = [0.25, 0.5, 0.25]
        SDDP.parameterize(sp, Ω, P) do ω
            p2 = ω + p1 * 0.8
            # Note the `-p2` becuase `fx2` appears on the right-hand side
            set_normalized_coefficient(c_cash, fx2, -p2)
            if t == 5
                set_objective_coefficient(sp, x2.out, p2)
            end
            return
        end
    end
    SDDP.train(model)
    return
end

demo()

thank you!This approach solved my problem.

[wrx1990]

Any other questions? If not, I will close this issue.

using SDDP, HiGHS, Test
function demo()
    model = SDDP.PolicyGraph(
    lattice,
    sense = :Max,
    upper_bound = 20.0,
    optimizer = HiGHS.Optimizer,
    ) do subproblem, node
        t, p1 = node::Tuple{Int,Float64}
        p2 = 1.0  # PLACEHOLDER
        @variable(sp, 0 <= x1, SDDP.State, initial_value = 0.0)
        @variable(sp, 0 <= x2, SDDP.State, initial_value = 0.0)
        @variable(sp, 0 <= x3, SDDP.State, initial_value = 10.0)
        @variables(sp, begin
            fx1     
            fx2     
            cash    
        end)
        @constraints(sp, begin
            x1.out == x1.in - fx1
            x2.out == x2.in - fx2
            c_cash, cash == p1 * fx1 + p2 * fx2
            x3.out == x3.in + cash 
        end)
        if t < 5
            @stageobjective(sp,  cash)
        else
            @stageobjective(sp, cash + p1 * x1.out + p2 * x2.out)
        end
        Ω = [-1.75, 0.08, 1.93]
        P = [0.25, 0.5, 0.25]
        SDDP.parameterize(sp, Ω, P) do ω
            p2 = ω + p1 * 0.8
            # Note the `-p2` becuase `fx2` appears on the right-hand side
            set_normalized_coefficient(c_cash, fx2, -p2)
            if t == 5
                set_objective_coefficient(sp, x2.out, p2)
            end
            return
        end
    end
    SDDP.train(model)
    return
end

demo()

In the code above, my variable p1 is obtained through a transition matrix. However, I want to modify the probabilities of the random variable ‘P’ to also be based on the transition probability matrix [[0.2, 0.4, 0.4], [0.3, 0.3, 0.4], [0.4, 0.1, 0.5]]. This means the current random value is dependent on the previous random value. How should I modify this part of the code?

[odow]

In the code above, my variable p1 is obtained through a transition matrix

How? It is the node in a lattice?

How should I modify this part of the code?

You can do stuff like:

if p1 > 0.5
     P = [0.25, 0.5, 0.25]
else
    P = [0.5, 0.5, 0.0]
end

[wrx1990]

How? It is the node in a lattice?

model = SDDP.PolicyGraph(
    lattice,
    sense = :Max,
    upper_bound = 20.0,
    optimizer = HiGHS.Optimizer,
    ) do subproblem, node
        t, p1 = node::Tuple{Int,Float64}
        p2 = 1.0  # PLACEHOLDER

The value of p1 in the above code is passed through lattice, and the results of lattice are as follows:

if p1 > 0.5
     P = [0.25, 0.5, 0.25]
else
    P = [0.5, 0.5, 0.0]
end```

I didn't express myself clearly, this part of the code doesn't meet my functional requirements.

Ω = [-1.75, 0.08, 1.93]
P = [0.25, 0.5, 0.25]
SDDP.parameterize(sp, Ω, P) do ω
      p2 = ω + p1 * 0.8

The current value of p2 is obtained through a random variable ω, with probabilities [0.25, 0.5, 0.25]. I now intend to change the probabilities of the random variable to transition probabilities, namely [[0.2, 0.4, 0.4], [0.3, 0.3, 0.4], [0.4, 0.1, 0.5]],as follows.

The probability of the random variable at time t is dependent on the value of the random variable at time t-1. How can I implement this functionality?

[wrx1990]

I am currently using the following method to solve this issue, not sure if it's correct.

using SDDP, HiGHS, Test
function demo()
    model = SDDP.PolicyGraph(
    lattice,
    sense = :Max,
    upper_bound = 20.0,
    optimizer = HiGHS.Optimizer,
    ) do subproblem, node
        t, p1 = node::Tuple{Int,Float64}
        p2 = 1.0  # PLACEHOLDER
        @variable(sp, 0 <= x1, SDDP.State, initial_value = 0.0)
        @variable(sp, 0 <= x2, SDDP.State, initial_value = 0.0)
        @variable(sp, 0 <= x3, SDDP.State, initial_value = 10.0)
        @variable(sp, random, SDDP.State, initial_value = 0.08)#Create a random state variable.

        @variables(sp, begin
            fx1     
            fx2     
            cash    
        end)
        @constraints(sp, begin
            x1.out == x1.in - fx1
            x2.out == x2.in - fx2
            c_cash, cash == p1 * fx1 + p2 * fx2
            x3.out == x3.in + cash 
        end)
        if t < 5
            @stageobjective(sp,  cash)
        else
            @stageobjective(sp, cash + p1 * x1.out + p2 * x2.out)
        end
        Ω = [-1.75, 0.08, 1.93]
        # Based on the previous random value, the current probability value is determined.
        if random.in==-1.75
            P = [0.2, 0.4, 0.4]
        elseif random.in==0.08
            P = [0.3, 0.3, 0.4]
        else
            P = [0.4, 0.1, 0.5]
        end
        
        SDDP.parameterize(sp, Ω, P) do ω
            p2 = ω + p1 * 0.8
            # Note the `-p2` becuase `fx2` appears on the right-hand side
            set_normalized_coefficient(c_cash, fx2, -p2)
            # Record random value
            JuMP.fix(random.out, ω)
            if t == 5
                set_objective_coefficient(sp, x2.out, p2)
            end
            return
        end
    end
    SDDP.train(model)
    return
end

demo()

[odow]

random.in==-1.75

You cannot do this.

Ignoring SDDP.jl, can you write out, in math, what your stochastic process is? How are p1 and p2 related?

You need to formulate your model as a policy graph, where the support and probability vectors depend only on the index of the node.

I don't fully understand your problem formulation, but you likely need to build a Markov chain with two dimensions: one for p1 and one for p2. If there are three possible values for p1 and three possible values for p2, then your Markov chain will have nine elements, so something like:

julia> P_I = [0.5 0.5 0.0; 0.5 0.0 0.5; 0.0 0.5 0.5]
3×3 Matrix{Float64}:
 0.5  0.5  0.0
 0.5  0.0  0.5
 0.0  0.5  0.5

julia> P_J = [0.2 0.4 0.4; 0.3 0.3 0.4; 0.4 0.1 0.5]
3×3 Matrix{Float64}:
 0.2  0.4  0.4
 0.3  0.3  0.4
 0.4  0.1  0.5

julia> indices = [(i, j) for i in 1:size(P_I, 1) for j in 1:size(P_J, 1)]
9-element Vector{Tuple{Int64, Int64}}:
 (1, 1)
 (1, 2)
 (1, 3)
 (2, 1)
 (2, 2)
 (2, 3)
 (3, 1)
 (3, 2)
 (3, 3)

julia> trannsition_matrix = [
           P1[i1, i2] * P2[j1, j2] for (i1, j1) in indices, (i2, j2) in indices
       ]
9×9 Matrix{Float64}:
 0.1   0.2   0.2   0.1   0.2   0.2   0.0   0.0   0.0
 0.15  0.15  0.2   0.15  0.15  0.2   0.0   0.0   0.0
 0.2   0.05  0.25  0.2   0.05  0.25  0.0   0.0   0.0
 0.1   0.2   0.2   0.0   0.0   0.0   0.1   0.2   0.2
 0.15  0.15  0.2   0.0   0.0   0.0   0.15  0.15  0.2
 0.2   0.05  0.25  0.0   0.0   0.0   0.2   0.05  0.25
 0.0   0.0   0.0   0.1   0.2   0.2   0.1   0.2   0.2
 0.0   0.0   0.0   0.15  0.15  0.2   0.15  0.15  0.2
 0.0   0.0   0.0   0.2   0.05  0.25  0.2   0.05  0.25

If there is correlation between I and J, then you'll need to account for that too.

[odow]

Any other questions?

[odow]

Closing because this seems resolved. Please comment if you have other questions.