Python gekko 方程定义中的换行符

我目前正在手动实现有限元的伽辽金法，并使用 python gekko 来求解所得的非线性代数方程组。这对于小型系统来说不会产生任何问题并且工作正常。

一旦系统变得更加复杂，涉及长方程和指数项m.exp()对于求解器来说，方程可能不再具有正确的格式。确切的错误消息是@error: Equation Definition Equation without an equality (=) or inequality (>,<).

我通过检查方程m.open_folder() in the model.apm。看来，以某种方式在方程中添加了换行符，从而失去了结果表达式之间的数学关系。这就是求解器无法识别的原因equality (=) or inequality (>,<)在第一个表达式中。 我还发现指数项m.exp(-EA / (R * T))导致断行。替换为np.exp(-EA/(R*T0))修复了该问题。这可能是由于 gekko 在后台缩短了长方程，并且此过程中存在错误。或者我的方程定义有问题。这是我的代码。

生成长而复杂的代数方程的基础

import pandas as pd
import numpy as np
from gekko import GEKKO
from typing import List, Callable
m = GEKKO(remote=False)
# constants
# -------------------------
R = 8.314  # J/mol K
# Rahmenbedingungen
# -------------------------
T0 = 550  # K
p0 = 15 * 10 ** 5 # Pa
# components
# -------------------------
x_bulk = np.array(
    [0.1, 0.08, 0.05]
)  # no unit Ethen, O2, CO2, EO, H2O
# catalyst pellet
# -------------------------
rho_p = 1171.2  # kg / m**3 Dichte eines Pellets
epsilon_p = 0.7  # Poroesitaet
tau = 4  # Tortuositaet
# heat transfer
# -------------------------
lambda_pm = 16  # W/m K Wärmeleitfähigkeit des reinen Katalysatorfestoffes
lambda_p = (1 - epsilon_p) * lambda_pm  # W/m K Wärmeleitfähigkeit des porösen Pellets
# reaction enthalpies
# -------------------------
dh_r1 = -107000  # J/mol
dh_r2 = -1323000  # J/mol
Deff = 5.42773893 * 10 **7

# reaction speeds 
# --------
def r_r1(p: np.array, T: float = T0) -> float:
    n_E = 0.6  # Reaktionsordnung Ethen
    n_O2 = 0.5  # Reaktionsordnung Sauerstoff
    k0 = 6.275 * 10**6  # pre-exponential collision factor [mol / (kg s Pa ** (nE+nO2))]
    EA = 74900 # activation energy [J/mol]
    K0 = 1.985 * 10**2  # pre-exponential adsorption factor [1/Pa]
    Tads = 2400  # adsorption temperature [K]
    r = k0 * m.exp(-EA / (R * T)) * p[0] ** n_E * p[1] ** n_O2 / (1 + K0 * np.exp(Tads / T0) * p[2])
    # this is m.exp() one cause for the quations not being formatted correctly
    # change it to np.exp(-EA/(R*T0)) and the equations are formatted correctly for the solver
    # this may be a mistake in my code or a bug in gekko
    return r

# differential equations 
# --------
def rhs(x: float, y: np.array, dy: np.array) -> np.array:
    # second derivatives of the partial pressures [bar/m**2]
    # use bar instead of Pa to increase numerical precision
    rhs_E   = -2 / x * dy[0] + R * y[-1] * T0 * rho_p / (Deff * p0) * (r_r1(p=y[0:-1] * p0, T=y[-1] * T0))
    rhs_O2  = -2 / x * dy[1] + R * y[-1] * T0 * rho_p / (Deff * p0) * (0.5 * r_r1(p=y[0:-1] * p0, T=y[-1] * T0))
    rhs_CO2 = -2 / x * dy[2] + R * y[-1] * T0 * rho_p / (Deff * p0) * (-2) * r_r1(p=y[0:-1] * p0, T=y[-1] * T0)
    rhs_T   = -2 / x * dy[-1] - rho_p / (lambda_p * T0) * (r_r1(p=y[0:-1] * p0, T=y[-1] * T0) * (-dh_r1))
    return np.array([rhs_E, rhs_O2, rhs_CO2, rhs_T])

# quadrature over the right hand side of the difffereantial equation
# ----------------

def give_quad_expressions(
    rhs: Callable, phi: List[np.poly1d], dphi: List[np.poly1d]
) -> Callable:
    """Creates an expression for the quadrature of the right-hand-side of the
    differential equation that is only dependent on the function values at the
    nodes y_i. It takes the quadrature base points to evaluate the
    base functions phi_i and their derivatives dphi_i to set up the expression.

    Args:
        rhs (Callable): Right hand side of the differential equation.
        Needs to be a function of x, y, and dy
        phi (List[np.poly1d]): List with all base polynomials phi. Mostly, base polynomials of Legendre type.
        dphi (List[np.poly1d]): List with the corresponding derivatives of the base polynomials.

    Returns:
        Callable: Expression for the weighted quadratures of the right-hand-side for each node that is only
        dependent on the unknown function values y_i. It can be called by a solving algorithm.
    """
    # these will be defined during runtime
    number_base_points = 2
    base_points = np.array([0.00021132, 0.00078868])
    weights = np.array([0.0005, 0.0005])

    # evaluate base functions at quadrature base points
    phi_eval = np.zeros(shape=(len(phi), number_base_points))
    dphi_eval = np.zeros(shape=(len(dphi), number_base_points))
    for index, phi_i in enumerate(phi):
        phi_eval[index] = phi_i(base_points)
        dphi_eval[index] = dphi[index](base_points)

    # set up the expression for the quadrature only dependent on the unknown y_i values
    def quad_expressions(y: np.array, model:GEKKO) -> np.array:
        """Variable expressions for the weighted quadratures of the integral of the weighted right-hand-side
        for each node that are only dependend on y.
        The integral is crucial for the Galerkin Method. int(phi_j * rhs(x, y, dy))dx.
        It allows the gekko solver to embed a numpy array y into the quadrature
        that is safed into the gekko model by using 'm.Array(m.Var)' for setting up the algebraic equation system.

        Args:
            y (np.array): Unknown function values that the gekko solver needs to find.

        Returns:
            np.array: Array with the weighted quadratures for each function and node as a function of y_i.
                first index specifies the function, second index the node
        """

        # TODO: Work with numpy arrays instead of lists to enhance speed.
        y_tilde = []
        dy_tilde = []
        # set up the trail functions and their derivatives using the base functions which are already
        # evaluated and the supposed y values at the nodes.
        for eq_l in range(y.shape[0]):
            y_tilde_eq = []
            dy_tilde_eq = []
            for bpoint_iq in range(number_base_points):
                # iterates over quadrature base points
                y_tilde_eq.append(np.dot(y[eq_l], phi_eval[:, bpoint_iq]))
                dy_tilde_eq.append(np.dot(y[eq_l], dphi_eval[:, bpoint_iq]))
            y_tilde.append(y_tilde_eq)
            dy_tilde.append(dy_tilde_eq)
        # evaluate the right-hand-sides with the base points and the supposed trial function values at the base points
        rhsides_evaluated_bpoints = rhs(
            base_points, np.array(y_tilde), np.array(dy_tilde)
        )
        # multiply the right-hand-side with the trial function of each node
        integrals = np.zeros_like(y)
        for eq_l, rhs_eq in enumerate(rhsides_evaluated_bpoints):
            # iterate over the equations
            for node_i in range(y.shape[1]):
                # iterate over the weight function of each node
                # calcualte the quadrature using the weights
                rhs_w_func = rhs_eq * phi_eval[node_i]
                dotp = model.sum(
                    [rhs_w_func[i] * weights[i] for i in range(number_base_points)]
                )
                # dotp = np.dot(
                #     rhs_w_func, self.weights
                # )
                integrals[eq_l, node_i] = dotp  # using the dot product via list comprehension with gekko to allow the model to split long equations
        return integrals

    return quad_expressions

建立方程

phi = [
    np.poly1d([ 2.e+06, -3.e+03,  1.e+00]),
    np.poly1d([-4000000.,     4000.,       -0.]),
    np.poly1d([ 2.e+06, -1.e+03,  0.e+00])
]

dphi = [
    np.poly1d([ 4.e+06, -3.e+03]),
    np.poly1d([-8.e+06,  4.e+03]),
    np.poly1d([ 4.e+06, -1.e+03])
]

diffusion_matrix = np.array(
    [[ 2333.33333333, -2666.66666667,   333.33333333],
    [-2666.66666667,  5333.33333333, -2666.66666667],
    [  333.33333333, -2666.66666667,  2333.33333333]]
)
y = m.Array(m.Var, (4, 3), value=x_bulk[0])

quad_expression = give_quad_expressions(rhs=rhs, phi=phi, dphi=dphi)
rhs_integrals = quad_expression(y=y, model=m)

for j in range(y.shape[0]):
    # iterates over each differential equation
    for i in range(y.shape[1]):
        # iterates over each node
        m.Equation(
            np.dot(
                y[j], (-1) * diffusion_matrix[:, i]
            ) == rhs_integrals[j, i]
        )

for eq in m._equations:
    print(eq)
m.solve()

Error

Exception: @error: Equation Definition

不带等式 (=) 或不等式 (>,(((((((v10)(-0.12200771519999987))+((v11)(0.6666554303999999)))+((v12)(0.4553522848000001))))(550)))))]))))([((((((((v1))(0.4553522848))+((v2)(0.6666554304000001)))+((v3)(-0.1220077152))))(1500000)))^(0.6)) 正在停止...

我发现了这个问题： 方程定义中出现这些换行符的原因是由于错误使用了 gekko 内置函数，例如m.exp(), m.sin()，...这些只能分段评估，而不是像我尝试的那样一次评估整个数组。

所以改变从

rhsides_evaluated_bpoints = rhs(
        base_points, np.array(y_tilde), np.array(dy_tilde)
    )

to

rhs_eval = []
    for i, bpoint in enumerate(base_points):
         rhs_eval.append(
              rhs(
                x=bpoint,
                y=y_tilde[:, i],
                dy=dy_tilde[:, i]
            )  
         )
    rhs_eval = np.array(rhs_eval).transpose()

完成了工作。