In [ ]:
import pfnet
import optalg
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [ ]:
T = 24
hours = range(T)
parser = pfnet.ParserMAT()
net = parser.parse('../../cases/case39.mat',T)
In [ ]:
for gen in net.generators:
gen.P_prev = gen.P_min
gen.dP_max = 0.3*(gen.P_max-gen.P_min)
In [ ]:
for load in net.loads:
load.P_max = 1.1*load.P
load.P_min = 0.9*load.P
load.target_power_factor = 0.9
In [ ]:
wind_profile = 0.5*(1+np.cos(np.linspace(0.,2.*np.pi,T)))
net.add_var_generators(net.get_load_buses(),50.,100.)
for gen in net.var_generators:
gen.P_ava = gen.P_ava[0]*wind_profile
gen.Q_max = 0.1*np.max(gen.P)
gen.Q_min = -0.1*np.max(gen.P)
plt.plot([sum(g.P_ava[t] for g in net.var_generators) for t in hours])
plt.ylabel('wind')
plt.axis('tight');
In [ ]:
net.add_batteries(net.get_generator_buses(),10.,40.)
In [ ]:
net.set_flags('bus',['variable','bounded'],'any','voltage magnitude')
net.set_flags('bus','variable','not slack','voltage angle')
net.set_flags('generator',['variable','bounded'],'any',['active power','reactive power'])
net.set_flags('variable generator',['variable','bounded'],'any',['active power','reactive power'])
net.set_flags('battery',['variable','bounded'],'any',['charging power','energy level'])
net.set_flags('load',['variable','bounded'],'any',['active power','reactive power'])
problem = pfnet.Problem(net)
problem.add_function(pfnet.Function('generation cost',1.,net))
problem.add_function(pfnet.Function('consumption utility',-1.,net))
problem.add_constraint(pfnet.Constraint('AC power balance',net))
problem.add_constraint(pfnet.Constraint('battery dynamics',net))
problem.add_constraint(pfnet.Constraint('variable bounds',net))
problem.add_constraint(pfnet.Constraint('load constant power factor',net))
problem.add_constraint(pfnet.Constraint('generator ramp limits',net))
problem.analyze()
problem.show()
In [ ]:
solver = optalg.opt_solver.OptSolverINLP()
solver.set_parameters({'quiet':True})
solver.solve(problem)
print(solver.get_status(), solver.get_iterations())
net.set_var_values(solver.get_primal_variables()[:net.num_vars])
In [ ]:
maxlp = np.max([sum(l.P[t] for l in net.loads) for t in hours])
lp = np.array([sum(l.P[t] for l in net.loads) for t in hours])/maxlp
gp = np.array([sum(g.P[t] for g in net.generators) for t in hours])/maxlp
ra = np.array([sum(g.P_ava[t] for g in net.var_generators) for t in hours])/maxlp
ru = np.array([sum(g.P[t] for g in net.var_generators) for t in hours])/maxlp
be = np.array([sum(b.E[t] for b in net.batteries) for t in hours])/maxlp
bc = np.array([sum(max([b.P[t],0.]) for b in net.batteries) for t in hours])/maxlp
bd = np.array([sum(max([-b.P[t],0.]) for b in net.batteries) for t in hours])/maxlp
# injections
fig1, ax1 = plt.subplots(1,1)
ax1.fill_between(hours,0,gp,facecolor='red')
ax1.fill_between(hours,gp,gp+bd,facecolor='magenta')
ax1.fill_between(hours,gp+bd,gp+bd+ru,facecolor='green')
ax1.fill_between(hours,gp+bd+ru,gp+bd+ra,facecolor='yellow')
ax1.plot([],[],color='red',label='generation')
ax1.plot([],[],color='magenta',label='battery')
ax1.plot([],[],color='green',label='renewables used')
ax1.plot([],[],color='yellow',label='renewabels left')
ax1.legend(loc='lower right')
ax1.set_ylabel('injection')
ax1.axis('tight')
# consumption
fig2, ax2 = plt.subplots(1,1)
ax2.fill_between(hours,0,lp,facecolor='blue')
ax2.fill_between(hours,lp,lp+bc,facecolor='cyan')
ax2.plot([],[],color='blue',label='load')
ax2.plot([],[],color='cyan',label='battery')
ax2.legend(loc='lower right')
ax2.set_ylabel('consumption')
ax2.axis('tight')
# load power factors
fig3, ax3 = plt.subplots(1,1)
for load in net.loads:
plt.plot(load.P,load.Q,'x')
ax3.set_xlabel('load active power')
ax3.set_ylabel('load reactive power')
ax3.axis('tight')
# battery energy
fig4, ax4 = plt.subplots(1,1)
for battery in net.batteries:
ax4.plot(hours,battery.E)
ax4.set_ylabel('battery energy')
ax4.axis('tight')
# battery charging
fig5, ax5 = plt.subplots(1,1)
for battery in net.batteries:
ax5.plot(hours,battery.P)
ax5.set_ylabel('battery charing power')
ax5.axis('tight');
In [ ]: