# <center><h1>Mathematical Formulation</h1></center>

The economical dispach of thermal generation units can be formulated as an objetive function, which should be minimized considering the limits of each power plant and also the supply of the demanded power. The following system of equations formulate the problem:

$$
\begin{equation}
    \begin{aligned}
       \text{Min } & F_t = \sum_{i=1}^{n_g} F_i(P_i) \\
       \text{s.t}  & \\
                   & \sum_{i=1}^{n_g} P_i = P_D \\
                   & P_i^{min} \leq P_i \leq P_i^{max} \\
    \end{aligned}
\end{equation}
$$

Onde:
 - $F_t$ = total financial cost of the dispach.
 - $F_i$ = financial cost function for the the generator $i$.
 - $P_i$ = power generated by the $i$. 
 - $P_D$ = power demand.
 - $P_i^{min}$ = minimum generation limit of the generator $i$.
 - $P_i^{max}$ = maximum generation limit of the generator $i$;
 - $n_g$ = number of generatio units.

A quadratic function can be used as a proxy of the fuel cost for each generation unit (Zhu, 2015). 

$$
\begin{equation}
    F_i(P_i) = c_iP_i^2 + b_iP_i + a_i
\end{equation}
$$

Bringin a more realistic approach to the economic dispach, a thermal generation unit might also have operative zones (Orero e Irving, 1996), which are directly related to auxiliary systems such as boilers and feed pumps. Such operative zones should be avoided as they bring risk to a safe operation of the plant, hence a discontinuous mathematical formulation rises:

$$
\begin{equation}
        P_i^{min} \leq P_i \leq P_i^{n_1} \\
        P_i^{n_2} \leq P_i \leq P_i^{n_3} \\
        P_i^{n_4} \leq P_i \leq P_i^{max} \\
\end{equation}
$$

---

# Solver

In order to solve the problem the [CVXOPT](https://cvxopt.org/) library will be used. The library offers a complete set of mathematical solvers implemented.

### Loading data

From now on we can use the `systems` module and initialize a `PowerSystem`, which should provide all of the necessary data for the solver. 


```python
from acopoweropt import system

PSystem = system.PowerSystem(name='s10')
```

### Solving

Defining a function to handle the solution of the PowerSystem


```python
from cvxopt import matrix, solvers

import pandas as pd
import numpy as np

def solve(demand:float, operation: pd.DataFrame,
                        maxiters: int=15,
                        show_progress: bool=False):
    """Returns a solution to a specific operation configuration

    Given a specific configuration to be solved, this function uses cvxopt
    quadratic programing to solve the system check:
        https://cvxopt.org/examples/tutorial/qp.html

    One possible source of configuration data can be obtained by using the
    provided `PowerSystem.sample_operation()` method, which randomly creates
    a possible configuration for the system to be operated

    Parameters
    ----------
    operation : pd.DataFrame
        DataFrame representing the operation of the system
    maxiters=15 : int
        Maximum number of iterations to be performed by the method.
    show_progress=False : bool
        Iteractively show progress during computation

    Returns
    -------
    dict
        A dictionary containing all of the solution results

    """

    # CVXOPT uses matrix like objects in order to model
    # a system of equations. Numpy can be used to prepare
    # the data so that the solver can be used.
    
    solvers.options['show_progress'] = show_progress
    # solvers.options['refinement'] = 2
    solvers.options['maxiters'] = maxiters

    demand = np.array([demand], dtype="double")

    # Equation parameters cP^2 + bP + a
    a = operation.a.sum()
    b = operation.b.to_numpy(dtype="double")
    c = operation.c.to_numpy(dtype="double")

    # CVXOPT needs a system of equations:
    Pmin = operation.Pmin.to_numpy(dtype="double")
    Pmax = operation.Pmax.to_numpy(dtype="double")

    P = matrix(2 * (c[..., None] * np.eye(operation.shape[0])))
    q = matrix(b)

    G_min = -1 * np.eye(operation.shape[0])
    G_max = np.eye(operation.shape[0])
    G = matrix(np.concatenate((G_min,G_max)))
    h = matrix(np.concatenate((-1 * Pmin,Pmax)))

    A = matrix(np.ones(operation.shape[0]), (1, c.shape[0]))
    b = matrix(demand)

    # Solving using Quadratic Programing
    solution = solvers.qp(P, q, G, h, A, b)

    # Interpreting the solution:
    Ft = solution.get('dual objective') + a
    Pg = pd.DataFrame([{'tgu': i+1,
                        'Pg': Pg} for i,Pg in enumerate(solution.get('x'))]
                      ).set_index('tgu')

    Fi = (((Pg.Pg ** 2) * operation.c) + (Pg.Pg * operation.b) + operation.c)
    Fi_df = pd.DataFrame({'tgu':Fi.index, 'Fi':Fi.values}).set_index('tgu')

    return {'status': solution.get('status'),
            'Ft': Ft,
            'operation': pd.concat([operation,Pg,Fi_df], axis=1)}
```

2ba26cc39316d96a24a7636b1239d25e44467642
ipynb
Jupyter Notebook
notebooks/SystemSolver.ipynb
ettoreaquino/acopoweropt
8c996358e17b92793c43eadb07db2ff77c5df516
MIT
2021-06-15T12:00:30.000Z
2021-06-15T12:00:30.000Z
notebooks/SystemSolver.ipynb
ettoreaquino/acopoweropt
8c996358e17b92793c43eadb07db2ff77c5df516
MIT
2021-06-14T01:58:41.000Z
2021-06-18T03:30:08.000Z
notebooks/SystemSolver.ipynb
ettoreaquino/acopoweropt
8c996358e17b92793c43eadb07db2ff77c5df516
MIT
Qwen/Qwen-72B
1. YES
2. YES