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

Topology of the power plant

Figure: Topology of the power plant

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_model and

  • to 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.

\[\eta_\mathrm{th}=\frac{|\sum P|}{\dot{Q}_{sg}}\]

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_objective method 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:

\[p_{e,1} > p_{e,2}\]

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
Scatter plot for all individuals during the optimization

Figure: Scatter plot for all individuals during the optimization

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()])