191 lines
5.7 KiB
Python
191 lines
5.7 KiB
Python
from pyomo.environ import *
|
|
from numpy import *
|
|
import matplotlib.pyplot as plt
|
|
|
|
# model
|
|
model = ConcreteModel()
|
|
|
|
# Set
|
|
model.t = Set(initialize=[1, 2, 3, 4], doc='time')
|
|
|
|
# Parameter
|
|
model.pb = Param(model.t, initialize={1: 3, 2: 1, 3: 2, 4: 4}, doc='Base load')
|
|
model.c = Param(model.t, initialize={1: 0.5, 2: 0.1, 3: 1, 4: 0.3}, doc='Cost ')
|
|
model.dt = Param(initialize=1, doc='time duration ')
|
|
# Variable
|
|
|
|
model.pl = Var(model.t,bounds=(0, None)) #todo example: 2
|
|
# model.pl = Var(RangeSet(1, 4),bounds=(0, 3)) #todo example: 3
|
|
# model.pl = Var(RangeSet(1, 4),bounds=(0, 3)) #todo example: 4
|
|
|
|
|
|
# model.pl = Var(RangeSet(1, 4),bounds=(0, 5)) #todo example: 5
|
|
# model.pmax = py.Param(initialize=5, doc='Maximum Power ')#todo example: 5
|
|
# model.emax = py.Param(initialize=8, doc='Maximum Controllable Energy ')#todo example: 5
|
|
# model.plmax = py.Param(initialize=6, doc='Maximum Controllable power')#todo example: 5
|
|
#
|
|
|
|
# Constraints
|
|
model.c1 = Constraint(expr=sum([model.pl[i] for i in model.t]) * model.dt == 8)
|
|
|
|
# todo this is for the example 5
|
|
# def total_power(model, i):
|
|
# return (model.pb[i] + model.pl[i]) <= model.pmax
|
|
# model.c2 = py.Constraint(model.t, rule=total_power, doc='Total power')
|
|
|
|
|
|
# Objective
|
|
|
|
model.obj = Objective(expr=sum([(model.pb[i] + model.pl[i]) * model.c[i] for i in model.t]) * model.dt)
|
|
|
|
# Solve
|
|
|
|
opt = SolverFactory('cplex')
|
|
result = opt.solve(model)
|
|
print(result)
|
|
# model.display()
|
|
# model.pprint()
|
|
|
|
# print variables
|
|
|
|
|
|
# print('------------Countable variable Solution------------------------')
|
|
for i in model.t:
|
|
print('x[%i] = %i' % (i, value(model.pl[i])))
|
|
print('objective function = ', value(model.obj))
|
|
|
|
# Plotting
|
|
time = []
|
|
pbase = []
|
|
pload = []
|
|
cost = []
|
|
|
|
for ll in range(len([model.t, model.pb, model.pl, model.c])):
|
|
if ll == 0:
|
|
for i in model.t:
|
|
time.append(value(model.t[i]))
|
|
elif ll == 1:
|
|
for i in model.pb:
|
|
pbase.append(value(model.pb[i]))
|
|
elif ll == 2:
|
|
for i in model.pl:
|
|
pload.append(value(model.pl[i]))
|
|
else:
|
|
for i in model.c:
|
|
cost.append(value(model.c[i]))
|
|
|
|
fig, (ax0, ax1) = plt.subplots(nrows=2,ncols= 1)
|
|
|
|
ax0.plot(time, pbase, 'b-', label='Base')
|
|
ax0.plot(time, pload, 'g--', label='Controable Load')
|
|
ptotal = [pbase[i] + pload[i] for i in range(len(pbase))]
|
|
ax0.plot(time, ptotal, 'r', label='Total Load')
|
|
ax0.set_ylabel('Electric power (W)')
|
|
ax0.legend(loc='upper left')
|
|
ax1.plot(time, cost, 'y-+', label='Cost')
|
|
ax1.legend(loc='upper left')
|
|
ax1.set_xlabel('Time (hour)')
|
|
ax1.set_ylabel('Electricity price ($/W/h)')
|
|
plt.suptitle("Energy management system-Example (4)")
|
|
plt.show()
|
|
|
|
|
|
'''
|
|
# another way do the coding, using excel sheet for the data and call them in the code
|
|
'''
|
|
|
|
|
|
# import matplotlib.pyplot as plt
|
|
# import pandas as pd
|
|
# import pyomo.environ as py
|
|
# # from pyomo.environ import *
|
|
#
|
|
# # model
|
|
# # from pyomo.core import value
|
|
#
|
|
# model = py.ConcreteModel()
|
|
# data = pd.read_excel('input_data.xlsx', sheet_name='set')
|
|
#
|
|
# # Set
|
|
# model.t = py.Set(initialize=data.time, doc='time')
|
|
#
|
|
# # Parameter
|
|
# model.pb = py.Param(model.t, initialize=dict(zip(data.time.values, data.base_load2.values)), doc='Base load')
|
|
# model.c = py.Param(model.t, initialize=dict(zip(data.time.values, data.price.values)), doc='cost')
|
|
# model.dt = py.Param(initialize=1, doc='time duration ')
|
|
# model.pmax = py.Param(initialize=5, doc='Maximum Power ')
|
|
# model.emax = py.Param(initialize=8, doc='Maximum Controllable Energy ')
|
|
# model.plmax = py.Param(initialize=6, doc='Maximum Controllable power')
|
|
# # for the simplicity of the equation
|
|
#
|
|
# # Variable
|
|
# model.pl = py.Var(model.t, bounds=(0, model.plmax))
|
|
#
|
|
#
|
|
#
|
|
# # Constraints
|
|
# # note:We can write the constraint like this
|
|
# # model.c1 = py.Constraint(expr=sum([pl[i] for i in t]) * dt == 8000)
|
|
#
|
|
# # note:Other way to write equations of constraints
|
|
# def Conctrol_power(model, i):
|
|
# return sum([model.pl[i] for i in model.t]) * model.dt == model.emax
|
|
# model.c1 = py.Constraint( rule=Conctrol_power, doc='Controllable power')
|
|
# def total_power(model, i):
|
|
# return (model.pb[i] + model.pl[i]) <= model.pmax
|
|
# model.c2 = py.Constraint(model.t, rule=total_power, doc='Total power')
|
|
#
|
|
#
|
|
# def objective_rule(model):
|
|
# return sum([(model.pb[i] + model.pl[i]) * model.c[i] for i in model.t])
|
|
# model.objective = py.Objective(rule=objective_rule, sense=py.minimize, doc='Define objective function')
|
|
#
|
|
# # Solve
|
|
# opt = py.SolverFactory('cplex')
|
|
# result = opt.solve(model)
|
|
# print(result)
|
|
# # model.pprint()
|
|
#
|
|
# # Print variables
|
|
# print('------------Controllable variable Solution------------------------')
|
|
# for i in model.t:
|
|
# print('x[%i] = %i' % (i, py.value(model.pl[i])))
|
|
# print('objective function = ', py.value(model.objective))
|
|
#
|
|
# # Plotting
|
|
#
|
|
# # Plotting
|
|
# time = []
|
|
# pbase = []
|
|
# pload = []
|
|
# cost = []
|
|
#
|
|
# for ll in range(len([model.t, model.pb, model.pl, model.c])):
|
|
# if ll == 0:
|
|
# for i in model.t:
|
|
# time.append(py.value(model.t[i]))
|
|
# elif ll == 1:
|
|
# for i in model.pb:
|
|
# pbase.append(py.value(model.pb[i]))
|
|
# elif ll == 2:
|
|
# for i in model.pl:
|
|
# pload.append(py.value(model.pl[i]))
|
|
# else:
|
|
# for i in model.c:
|
|
# cost.append(py.value(model.c[i]))
|
|
#
|
|
# fig, (ax0, ax1) = plt.subplots(nrows=2, ncols=1)
|
|
#
|
|
# ax0.plot(time, pbase, 'b-', label='Base')
|
|
# ax0.plot(time, pload, 'g--', label='Controable Load')
|
|
# ptotal = [pbase[i] + pload[i] for i in range(len(pbase))]
|
|
# ax0.plot(time, ptotal, 'r', label='Total Load')
|
|
# ax0.set_ylabel('Electric power (W)')
|
|
# ax0.legend(loc='upper left')
|
|
# ax1.plot(time, cost, 'y-+', label='Cost')
|
|
# ax1.legend(loc='upper left')
|
|
# ax1.set_xlabel('Time (hour)')
|
|
# ax1.set_ylabel('Electricity price ($/W/h)')
|
|
# plt.suptitle("Energy management system-Example (5)")
|
|
# plt.show()
|