ProblemFormulation.py

# load the pyomo library
from pyomo.core import *

# model and network definition
model = AbstractModel()
model.Nodes=Set(ordered=True)
model.Lines=Set(within=model.Nodes*model.Nodes)
model.Types=Param(model.Nodes,default=0) #2 PV, 0 PQ, 3 Slack

## definition of parameters, a high default is a safe choice to avoid unwanted behaviours
# power flow parameters
model.X=Param(model.Lines,default=1000000)
model.Costs=Param(model.Nodes,default=0)
model.P_gen=Param(model.Nodes)
model.P_load=Param(model.Nodes)
model.P_gen_min=Param(model.Nodes)
model.P_gen_max=Param(model.Nodes)

# problem-related parameters
model.bmva=Param()
model.ccv=Param(model.Lines,default=1000000)
model.Plines_max1=Param(model.Lines)
def P_lines_max_initialize(model,n,m):
	return model.Plines_max1[n,m]/1
model.Plines_max=Param(model.Lines,initialize=P_lines_max_initialize)
#stochastic parameters
model.cpvB=Param(model.Nodes,default=1000000)
model.cpvC=Param(model.Nodes,default=1000000)
model.P_fluctB=Param(model.Nodes,default=0)
model.P_fluctC=Param(model.Nodes,default=0)

#scenarios' variables
def a0_initialize(model,n,m):
	if n=='W' and m=='o' :
		return 1
	elif n=='M' and m=='m':
		return 1
	elif n=='G' and m=='a':
		return 1
	else :
		model.a0[n,m].fixed=True
	return 1	
	
#model.a0=Var(model.Lines,bounds=(0.5,1),initialize=a0_initialize)
#model.a0=Param(model.Lines,default=1)
model.a0=Var(model.Lines,within=Binary,initialize=a0_initialize) #this works only with cbc

def init_P_deltaB(model,n):
	if model.Types[n]==0 and model.P_gen[n]==0 :
		model.P_deltaB[n].fixed=True
	return 0
def init_P_deltaC(model,n):
	if model.Types[n]==0 and model.P_gen[n]==0 :
		model.P_deltaC[n].fixed=True
	return 0

model.P_deltaB=Var(model.Nodes,initialize=init_P_deltaB,within=NonNegativeReals)#,bounds=init_bounds_P_deltaB)
model.P_deltaC=Var(model.Nodes,initialize=init_P_deltaC,within=NonNegativeReals)#,bounds=init_bounds_P_deltaC)

#auxiliary variables
def init_anglesA(model,n):
	if model.Types[n]==3 :
		model.anglesA[n].fixed=True
	return 0
def init_anglesB(model,n):
	if model.Types[n]==3 :
		model.anglesB[n].fixed=True
	return 0
def init_anglesC(model,n):
	if model.Types[n]==3 :
		model.anglesC[n].fixed=True
	return 0
model.anglesA=Var(model.Nodes,initialize=init_anglesA,bounds=(-3.14,3.14))
model.anglesB=Var(model.Nodes,initialize=init_anglesB,bounds=(-3.14,3.14))
model.anglesC=Var(model.Nodes,initialize=init_anglesC,bounds=(-3.14,3.14))

#power flow constraints
def powerFlow_firstStage_rule(model,i):
	totalFlow = model.P_gen[i]-model.P_load[i]
	totalFlow += sum((model.anglesA[n]-model.anglesA[m])/model.X[n,m] for (n,m) in model.Lines if n==i and model.a0[n,m]!=0)
	totalFlow -= sum((model.anglesA[n]-model.anglesA[m])/model.X[n,m] for (n,m) in model.Lines if m==i and model.a0[n,m]!=0)
	return (-0.01,totalFlow,0.01) 
	#return totalFlow==0 uncomment if more accuracy is required, it may not converge
def powerFlow_secondStage_rule(model,i):
	totalFlow = model.P_gen[i]+model.P_fluctB[i]/model.bmva-model.P_deltaB[i]-model.P_load[i]
	totalFlow += sum((model.anglesB[n]-model.anglesB[m])/model.X[n,m] for (n,m) in model.Lines if n==i and model.a0[n,m]!=0)
	totalFlow -= sum((model.anglesB[n]-model.anglesB[m])/model.X[n,m] for (n,m) in model.Lines if m==i and model.a0[n,m]!=0)
	return (-0.01,totalFlow,0.01)
	#return totalFlow==0
def powerFlow_thirdStage_rule(model,i):
	totalFlow = model.P_gen[i]+model.P_fluctC[i]/model.bmva-model.P_deltaC[i]-model.P_load[i]
	totalFlow += sum((model.anglesC[n]-model.anglesC[m])/model.X[n,m] for (n,m) in model.Lines if n==i and model.a0[n,m]!=0)
	totalFlow -= sum((model.anglesC[n]-model.anglesC[m])/model.X[n,m] for (n,m) in model.Lines if m==i and model.a0[n,m]!=0)
	return (-0.01,totalFlow,0.01)
	#return totalFlow==0
		
model.PowerFlowConstraintA = Constraint(model.Nodes,rule=powerFlow_firstStage_rule)
model.PowerFlowConstraintB = Constraint(model.Nodes,rule=powerFlow_secondStage_rule)
model.PowerFlowConstraintC = Constraint(model.Nodes,rule=powerFlow_thirdStage_rule)

#max line flow constraint
def max_line_flowA_rule(model,n,m):
	return (-1,(model.anglesA[n]-model.anglesA[m])*model.a0[n,m]/model.X[n,m]/(model.Plines_max[n,m]/model.bmva),1)
def max_line_flowB_rule(model,n,m):
	return (-1,(model.anglesB[n]-model.anglesB[m])*model.a0[n,m]/model.X[n,m]/(model.Plines_max[n,m]/model.bmva),1)
def max_line_flowC_rule(model,n,m):
	return (-1,(model.anglesC[n]-model.anglesC[m])*model.a0[n,m]/model.X[n,m]/(model.Plines_max[n,m]/model.bmva),1)

model.MaxLineFlowMaxA = Constraint(model.Lines,rule=max_line_flowA_rule)
model.MaxLineFlowMaxB = Constraint(model.Lines,rule=max_line_flowB_rule)
model.MaxLineFlowMaxC = Constraint(model.Lines,rule=max_line_flowC_rule)


#max deltaP limit
def max_DeltaP_limit(model,n):
	if model.Types[n]!=0:
		return (0,model.P_deltaB[n]+model.P_deltaC[n]-(model.P_fluctB[n]+model.P_fluctC[n])/model.bmva,model.P_gen[n])
	else :
		return Constraint.Skip
model.maxDeltaP = Constraint(model.Nodes,rule=max_DeltaP_limit)


#Stage-specific cost computations
def ComputeFirstStageCost_rule(model):
    return sum((1-model.a0[n,m])*model.ccv[n,m] for (n,m) in model.Lines)
def ComputeSecondStageCost_rule(model):
    return sum(model.cpvB[n]*model.P_deltaB[n] for n in model.Nodes if model.Types[n]!=0)
def ComputeThirdStageCost_rule(model):
	return sum(model.cpvC[n]*model.P_deltaC[n] for n in model.Nodes if model.Types[n]!=0)

model.FirstStageCost = Expression(rule=ComputeFirstStageCost_rule)
model.SecondStageCost = Expression(rule=ComputeSecondStageCost_rule)
model.ThirdStageCost = Expression(rule=ComputeThirdStageCost_rule)


# obj. function -> minimize: sum of StageCosts

def total_cost_rule(model):
    return model.FirstStageCost + model.SecondStageCost + model.ThirdStageCost

model.Total_Cost_Objective = Objective(rule=total_cost_rule, sense=minimize)