Thermal Power Plant Efficiency Optimization#
Task#
Designing a power plant meets multiple different tasks, such as finding the optimal fresh steam temperature and pressure to reduce exhaust steam water content, or the optimization of extraction pressures to maximize cycle efficiency and many more.
In case of a rather simple power plant topologies the task of finding optimized values for e.g. extraction pressures is still manageable without any optimization tool. As the topology becomes more complex and boundary conditions come into play the usage of additional tools is recommended. The following tutorial is intended to show the usage of PyGMO in combination with TESPy to maximize the cycle efficiency of a power plant with two extractions.
You can download the code here:
optimization_example.py
Figure: Topology of the power plant#
Figure: Topology of the power plant#
What is PyGMO?#
PyGMO (Python Parallel Global Multiobjective Optimizer, [3]) is a library that provides numerous evolutionary optimization algorithms. PyGMO can be used to solve constrained, unconstrained, single objective and multi objective problems.
Evolutionary Algorithms#
Evolutionary Algorithms (EA) are optimization algorithms inspired by biological evolution. In a given population the algorithm uses the so-called fitness function to determine the quality of the solutions to each individual (set of decision variables) problem. The best possible solution of the population is called champion. Via mutation, recombination and selection your population evolves to find better solutions.
EA will never find an exact solution to your problem. They can only give an approximation for the real optimum.
Install PyGMO#
With the conda package manager PyGMO is available for Linux, OSX and Windows thanks to the infrastructure of conda-forge:
conda install -c conda-forge pygmo
Windows user can perform an installation from source as an alternative to conda. For further information on this process we recommend the PyGMO installation accordingly.
On Linux you also have the option to use the pip package installer:
pip install pygmo
Creating your TESPy-Model#
To use the API, you need to define a class holding a TESPy network. The initialization of the class should build the network and run an initial simulation. Furthermore, you have to define methods
to get component or connection parameters of the plant
get_param,to run a new simulation for every new input from PyGMO
solve_modelandto return the objective value
get_objective.
First, we set up the class with the TESPy network.
Display source code for the PowerPlant class
import numpy as np
import pygmo as pg
from tespy.components import CycleCloser
from tespy.components import Sink
from tespy.components import Source
from tespy.components import Condenser
from tespy.components import Desuperheater
from tespy.components import SimpleHeatExchanger
from tespy.components import Merge
from tespy.components import Splitter
from tespy.components import Pump
from tespy.components import Turbine
from tespy.connections import Bus
from tespy.connections import Connection
from tespy.networks import Network
from tespy.tools.optimization import OptimizationProblem
class SamplePlant:
"""Class template for TESPy model usage in optimization module."""
def __init__(self):
self.nw = Network()
self.nw.set_attr(
p_unit="bar", T_unit="C", h_unit="kJ / kg", iterinfo=False
)
# components
# main cycle
sg = SimpleHeatExchanger("steam generator")
cc = CycleCloser("cycle closer")
hpt = Turbine("high pressure turbine")
sp1 = Splitter("splitter 1", num_out=2)
mpt = Turbine("mid pressure turbine")
sp2 = Splitter("splitter 2", num_out=2)
lpt = Turbine("low pressure turbine")
con = Condenser("condenser")
pu1 = Pump("feed water pump")
fwh1 = Condenser("feed water preheater 1")
fwh2 = Condenser("feed water preheater 2")
dsh = Desuperheater("desuperheater")
me2 = Merge("merge2", num_in=2)
pu2 = Pump("feed water pump 2")
pu3 = Pump("feed water pump 3")
me = Merge("merge", num_in=2)
# cooling water
cwi = Source("cooling water source")
cwo = Sink("cooling water sink")
# connections
# main cycle
c0 = Connection(sg, "out1", cc, "in1", label="0")
c1 = Connection(cc, "out1", hpt, "in1", label="1")
c2 = Connection(hpt, "out1", sp1, "in1", label="2")
c3 = Connection(sp1, "out1", mpt, "in1", label="3", state="g")
c4 = Connection(mpt, "out1", sp2, "in1", label="4")
c5 = Connection(sp2, "out1", lpt, "in1", label="5")
c6 = Connection(lpt, "out1", con, "in1", label="6")
c7 = Connection(con, "out1", pu1, "in1", label="7", state="l")
c8 = Connection(pu1, "out1", fwh1, "in2", label="8", state="l")
c9 = Connection(fwh1, "out2", me, "in1", label="9", state="l")
c10 = Connection(me, "out1", fwh2, "in2", label="10", state="l")
c11 = Connection(fwh2, "out2", dsh, "in2", label="11", state="l")
c12 = Connection(dsh, "out2", me2, "in1", label="12", state="l")
c13 = Connection(me2, "out1", sg, "in1", label="13", state="l")
self.nw.add_conns(
c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13
)
# preheating
c21 = Connection(sp1, "out2", dsh, "in1", label="21")
c22 = Connection(dsh, "out1", fwh2, "in1", label="22")
c23 = Connection(fwh2, "out1", pu2, "in1", label="23")
c24 = Connection(pu2, "out1", me2, "in2", label="24")
c31 = Connection(sp2, "out2", fwh1, "in1", label="31")
c32 = Connection(fwh1, "out1", pu3, "in1", label="32")
c33 = Connection(pu3, "out1", me, "in2", label="33")
self.nw.add_conns(c21, c22, c23, c24, c31, c32, c33)
# cooling water
c41 = Connection(cwi, "out1", con, "in2", label="41")
c42 = Connection(con, "out2", cwo, "in1", label="42")
self.nw.add_conns(c41, c42)
# busses
# power bus
self.power = Bus("power")
self.power.add_comps(
{"comp": hpt, "char": -1}, {"comp": mpt, "char": -1},
{"comp": lpt, "char": -1}, {"comp": pu1, "char": -1},
{"comp": pu2, "char": -1}, {"comp": pu3, "char": -1}
)
# heating bus
self.heat = Bus("heat")
self.heat.add_comps({"comp": sg, "char": 1})
self.nw.add_busses(self.power, self.heat)
hpt.set_attr(eta_s=0.9)
mpt.set_attr(eta_s=0.9)
lpt.set_attr(eta_s=0.9)
pu1.set_attr(eta_s=0.8)
pu2.set_attr(eta_s=0.8)
pu3.set_attr(eta_s=0.8)
sg.set_attr(pr=0.92)
con.set_attr(pr1=1, pr2=0.99, ttd_u=5)
fwh1.set_attr(pr1=1, pr2=0.99, ttd_u=5)
fwh2.set_attr(pr1=1, pr2=0.99, ttd_u=5)
dsh.set_attr(pr1=0.99, pr2=0.99)
c1.set_attr(m=200, T=650, p=100, fluid={"water": 1})
c2.set_attr(p=20)
c4.set_attr(p=3)
c41.set_attr(T=20, p=3, fluid={"INCOMP::Water": 1})
c42.set_attr(T=28, p0=3, h0=100)
# parametrization
# components
self.nw.solve("design")
self.stable = "_stable"
self.nw.save(self.stable)
self.solved = True
self.nw.print_results()
Next, we add the methods get_param, solve_model and
get_objective. On top of that, we add a setter working similarly as the
getter. The objective is to maximize thermal efficiency as defined in the
equation below.
Attention
The sense of optimization is always minimization, therefore you need to define your objective functions in the appropriate way: We return the reciprocal value of the efficiency for this reason.
We also have to make sure, only the results of physically feasible solutions
are returned. In case we have infeasible solutions, we can simply return
np.nan. An infeasible solution is obtained in case the power of a
turbine is positive, the power of a pump is negative, or the heat exchanged
in any of the preheaters is positive. We also check, if the calculation does
converge.
Display source code for the class methods
def get_param(self, obj, label, parameter):
"""Get the value of a parameter in the network"s unit system.
Parameters
----------
obj : str
Object to get parameter for (Components/Connections).
label : str
Label of the object in the TESPy model.
parameter : str
Name of the parameter of the object.
Returns
-------
value : float
Value of the parameter.
"""
if obj == "Components":
return self.nw.get_comp(label).get_attr(parameter).val
elif obj == "Connections":
return self.nw.get_conn(label).get_attr(parameter).val
def set_params(self, **kwargs):
if "Connections" in kwargs:
for c, params in kwargs["Connections"].items():
self.nw.get_conn(c).set_attr(**params)
if "Components" in kwargs:
for c, params in kwargs["Components"].items():
self.nw.get_comp(c).set_attr(**params)
def solve_model(self, **kwargs):
"""
Solve the TESPy model given the the input parameters
"""
self.set_params(**kwargs)
self.solved = False
try:
self.nw.solve("design")
if not self.nw.converged:
self.nw.solve("design", init_only=True, init_path=self.stable)
else:
# might need more checks here!
if (
any(self.nw.results["Condenser"]["Q"] > 0)
or any(self.nw.results["Desuperheater"]["Q"] > 0)
or any(self.nw.results["Turbine"]["P"] > 0)
or any(self.nw.results["Pump"]["P"] < 0)
):
self.solved = False
else:
self.solved = True
except ValueError as e:
self.nw.lin_dep = True
self.nw.solve("design", init_only=True, init_path=self.stable)
def get_objective(self, objective=None):
"""
Get the current objective function evaluation.
Parameters
----------
objective : str
Name of the objective function.
Returns
-------
objective_value : float
Evaluation of the objective function.
"""
if self.solved:
if objective == "efficiency":
return 1 / (
self.nw.busses["power"].P.val /
self.nw.busses["heat"].P.val
)
else:
msg = f"Objective {objective} not implemented."
raise NotImplementedError(msg)
else:
return np.nan
After this, we import the
tespy.tools.optimization.OptimizationProblem class and create an
instance of our self defined class, which we pass to an instance of the
OptimizationProblem class. We also have to pass
the variables to optimize,
the constraints to consider and
the objective function name (you could define multiple in the
get_objectivemethod if you wanted).
We set one inequality constraint, namely that the pressure of the first extraction has to be higher than the pressure at the second one:
To do this, we can set a lower limit for the pressure at connection 2 and reference the pressure at connection 4 as seen in the code:
plant = SamplePlant()
plant.get_objective("efficiency")
variables = {
"Connections": {
"2": {"p": {"min": 1, "max": 40}},
"4": {"p": {"min": 1, "max": 40}}
}
}
constraints = {
"lower limits": {
"Connections": {
"2": {"p": "ref1"}
},
},
"ref1": ["Connections", "4", "p"]
}
optimize = OptimizationProblem(
plant, variables, constraints, objective="efficiency"
)
Before we can run the optimization, we only need to select an appropriate
algorithm. After that we can start the optimization run. For more information
on algorithms available in the PyGMO framework and their individual
specifications please refer to the respective section in their online
documentation:
list of algorithms.
Create an initial population and then specify the number of individuals, the
number of generations and call the
tespy.tools.optimization.OptimizationProblem.run() method of your
OptimizationProblem instance passing the algorithm, the population and
the number of individuals and generations.
Run PyGMO-Optimization#
The following code then simply runs the PyGMO optimization.
num_ind = 10
num_gen = 100
# for algorithm selection and parametrization please consider the pygmo
# documentation! The number of generations indicated in the algorithm is
# the number of evolutions we undertake within each generation defined in
# num_gen
algo = pg.algorithm(pg.ihs(gen=3, seed=42))
# create starting population
pop = pg.population(pg.problem(optimize), size=num_ind, seed=42)
optimize.run(algo, pop, num_ind, num_gen)
In our run, we got:
Efficiency: 44.82 %
Extraction 1: 25.4756 bar
Extraction 2: 2.5428 bar
Figure: Scatter plot for all individuals during the optimization#
Figure: Scatter plot for all individuals during the optimization#
Finally, you can access the individuals in each of the generations, and you can have a look at you population. For more info on the population API please visit the pygmo documentation.
# To access the results
print(optimize.individuals)
# check pygmo documentation to see, what you can get from the population
pop
# plot the results
import matplotlib.pyplot as plt
# make text reasonably sized
plt.rc("font", **{"size": 18})
fig, ax = plt.subplots(1, figsize=(16, 8))
filter_valid_constraint = optimize.individuals["valid"].values
filter_valid_result = ~np.isnan(optimize.individuals["efficiency"].values)
data = optimize.individuals.loc[filter_valid_constraint & filter_valid_result]
sc = ax.scatter(
data["Connections-2-p"],
data["Connections-4-p"],
c=1 / data["efficiency"] * 100,
s=100
)
cbar = plt.colorbar(sc)
cbar.set_label("Thermal efficiency in %")
ax.set_axisbelow(True)
ax.set_xlabel("Pressure at connection 2 in bar")
ax.set_ylabel("Pressure at connection 4 in bar")
plt.tight_layout()
fig.savefig("pygmo_optimization.svg")
print(data.loc[data["efficiency"].values == data["efficiency"].min()])