Strong Duality constraint using Pyomo library in Python

Question

I am implementing an Inverse Optimization Problem which instead of using KKT conditions applies the strong duality theorem (primal optimal objective and the dual optimal objective are equal). For this purpose, it is required to formulate the primal and the dual problem. I have 4 generators with 8 energy blocks (b) each and a fixed demand. This is a market clearing per hour. The primal problem is an easy Market Clearing and it is as follows using the library Pyomo:

model = ConcreteModel()
model.g1=Var(b, within=NonNegativeReals)
model.g2=Var(b, within=NonNegativeReals)
model.g3=Var(b, within=NonNegativeReals)
model.g4=Var(b, within=NonNegativeReals)

model.obj = Objective(expr=
                      (sum(g1price[i]*model.g1[i] for i in b)+
                       sum(g2price[i]*model.g2[i] for i in b)+
                       sum(g3price[i]*model.g3[i] for i in b)+
                       sum(g4price[i]*model.g4[i] for i in b)))
model.con_power_balance=Constraint(
                      expr=
                      (sum(model.g1[i] for i in b)+
                       sum(model.g2[i] for i in b)+
                       sum(model.g3[i] for i in b)+
                       sum(model.g4[i] for i in b)- demand) == 0)


model.con_g1max=ConstraintList()
for i in b:
    model.con_g1max.add(model.g1[i] <= gsize[i])
model.con_g2max=ConstraintList()
for i in b:
    model.con_g2max.add(model.g2[i] <= gsize[i])
model.con_g3max=ConstraintList()
for i in b:
    model.con_g3max.add(model.g3[i]< = gsize[i])
model.con_g4max=ConstraintList()
for i in b:
    model.con_g4max.add(model.g4[i] <= gsize[i])

Let's say that the previous equations are named (1a-1f). When no stated, is set as default to minimize the objective function. As the primal problem is a minimization problem, its dual must be a maximization one. But it is equivalent max[F(x)] == min[-F(x)], and it is what I have applied. Here is the dual problem (equations 2a-2e):

model = ConcreteModel()
model.mu_g1max=Var(b, within=NonNegativeReals)
model.mu_g2max=Var(b, within=NonNegativeReals)
model.mu_g3max=Var(b, within=NonNegativeReals)
model.mu_g4max=Var(b, within=NonNegativeReals)
model.mu_g1min=Var(b, within=NonNegativeReals)
model.mu_g2min=Var(b, within=NonNegativeReals)
model.mu_g3min=Var(b, within=NonNegativeReals)
model.mu_g4min=Var(b, within=NonNegativeReals)
model.lambda=Var(b, within=NonNegativeReals)

model.obj = Objective(expr=
                      (sum(gsize[i]*model.mu_g1max[i] for i in b)+
                       sum(gsize[i]*model.mu_g2max[i] for i in b)+
                       sum(gsize[i]*model.mu_g3max[i] for i in b)+
                       sum(gsize[i]*model.mu_g4max[i] for i in b))

model.con_g1_dual=ConstraintList()
for i in b:
    model.con_g1_dual.add(model.lambda[i]+model.mu_g1min[i]-model.mu_g1max <= gsize[i])
model.con_g2_dual=ConstraintList()
for i in b:
    model.con_g2_dual.add(model.lambda[i]+model.mu_g2min[i]-model.mu_g2max <= gsize[i])
model.con_g3_dual=ConstraintList()
for i in b:
    model.con_g3_dual.add(model.lambda[i]+model.mu_g3min[i]-model.mu_g3max <= gsize[i])
model.con_g4_dual=ConstraintList()
for i in b:
    model.con_g4_dual.add(model.lambda[i]+model.mu_g4min[i]-model.mu_g4max <= gsize[i])

Once the primal and the dual problem are presented, the real problem that I have to solve is the following one:

model = ConcreteModel()
model.lambda=Var(b, within=NonNegativeReals)
model.obj = Objective(expr=
                      np.absolute(sum(model.lambda[i]-lambda_ini[i] for i in b*4)))

Here I have to put as constraints the ones from the dual problem (2b-2e) written before and moreover apply the strong duality constraint which would be:

min (objective function primal problem) = max (objective function dual problem)

Here is my problem: How would be this last constraint written??

What you are describing ought not to be very difficult. For one, you can obtain dual values using the model.dual = Suffix(direction=Suffix.IMPORT_EXPORT) suffix. I would not use a variable named lambda, simply because lambda in programming can have a very different meaning. — Qi Chen
– Qi Chen, Commented Dec 5, 2017 at 15:24
I am a bit confused what you mean when you refer to equations 2a-2e and 2b-2e, because they are not labeled in your Pyomo formulation as thus. — Qi Chen
– Qi Chen, Commented Dec 5, 2017 at 15:25

oakca · Accepted Answer · 2017-12-21 14:56:35Z

1

Well, I did not really understand your question, and I cannot comment yet, so I am posting an answer, hoping that it would help you.

First of all in pyomo, like @Qi Chen mentioned, you can access the dual values which generated by solver itself. For that you need to add following in your model:

model.dual = Suffix(direction=Suffix.IMPORT)

Then you may call the instance of your model elsewhere, and you can get the dual values of the constraints, which are inside your model via:

dual_value = inst.dual[inst.constraint_name]

e.g.:

dual_con_g1max = inst.dual[inst.con_g1max]

but, ofc since con_g1max is a Constraint_List, you need to give its index to the dual getting algo. like:

inst.dual[inst.con_g1max[1]]

answered Dec 21, 2017 at 14:56

oakca

1,5883 gold badges19 silver badges48 bronze badges

Sign up to request clarification or add additional context in comments.

Collectives™ on Stack Overflow

Strong Duality constraint using Pyomo library in Python

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related