Python和MATLAB锂电铅蓄电化学微分模型和等效电路

🎯要点

1. 对比三种电化学颗粒模型：电化学的锂离子电池模型、单粒子模型和带电解质的单粒子模型。
2. 求解粒子域内边界通量与局部电流密度有关的扩散方程。
3. 扩展为两个相的负或正电极复合电极粒子模型。
4. 模拟四种耦合机制下活性物质损失情况。
5. 模拟锂离子电池三参数模型等效电路。

🍁电池解析模型

🍪语言内容分比

🍇Python和MATLAB符号微分

计算导数的常用方法有 3 种：

• 数值微分
• 符号微分
• 自动微分

数值微分依赖于导数的定义：
$$f^{\prime }(x)=\underset{h\to 0}{\lim }\frac{f(x+h)-f(x)}{h}$$
其中，你把一个非常小的 $h$ 和评估函数放在两个地方。这是最基本的公式，在实践中，人们使用其他公式，这些公式给出的估计误差更小。这种计算导数的方法主要适用于你不知道你的函数并且只能对其进行采样的情况。此外，对于高维函数，它需要大量的计算。

符号微分操纵数学表达式。如果你曾经使用过 matlab 或 mathematica，那么你会看到类似这样的内容

$$\begin{array}{rl}& \frac{d}{dx}\left(\frac{x^2\cos (x-7)}{\sin (x)}\right)=\\ & x^2\sin (7-x)\csc (x)+x^2(-\cos (7-x))\cot (x)\csc (x)+2x\cos (7-x)\csc (x)\end{array}$$
在这里，他们知道每个数学表达式的导数，并使用各种规则（乘积规则、链式规则）来计算最终的导数。然后他们简化表达式以获得最终表达式。

自动微分可操纵计算机程序块。微分器具有对程序的每个元素求导的规则。它还使用链式法则将复杂表达式分解为更简单的表达式。

Python符号微分

导数规则要求一个项有 2 个子项，即二元运算（例如 *），以及一个一元运算（例如 sin）。

term.py:

from functools import reduce

class Term:

    def __init__(self, label: str or int, children=None):
        self.label = label
        self.children = children if children else []

    def __repr__(self):
        if self.children:
            return f"({self.label} " + reduce(lambda str1, str2: str1+" "+str2, (repr(child) for child in self.children)) + ")"
        else:
            return f"{self.label}"

    @classmethod
    def copy(cls, obj):
        """makes a copy of obj"""
        return cls(obj.label, [cls.copy(child) for child in obj.children])

symbols_table:

MUL = '*'
ADD = '+'
SUB = '-'
DIV = '/'
POW = '^'
LOG = 'log'
SIN = 'sin'
COS = 'cos'
TAN = 'tan'
OPEN_BR = '('
CLOSE_BR = ')'
VAR = 'x'

symbols = {LOG, SIN, COS, TAN, MUL, ADD, SUB, DIV, POW, OPEN_BR, CLOSE_BR, VAR}

from term import Term
from symbols_table import *
import numbers

def calculate_derivative(term1: Term, rules: dict) -> Term:
    if isinstance(term1.label, numbers.Number):
        return Term(0)
    elif term1.label == VAR:
        return Term(1)
    else:
        return rules[term1.label](term1, rules)


def mul_derivative(term1: Term, rules) -> Term:
    child1, child2 = term1.children  
    return Term(ADD, [Term(MUL, [calculate_derivative(child1, rules), child2]), Term(MUL, [child1, calculate_derivative(child2, rules)])])


def div_derivative(term1: Term, rules) -> Term:
    child1, child2 = term1.children  
    return Term(DIV, [Term(SUB, [Term(MUL, [calculate_derivative(child1, rules), child2]), Term(MUL, [child1, calculate_derivative(child2, rules)])]),
                      Term(POW, [child2, Term(2)])])


def sum_derivative(term1: Term, rules) -> Term:
    child1, child2 = term1.children  
    return Term(term1.label, [calculate_derivative(child1, rules), calculate_derivative(child2, rules)])


def power_derivative(term1: Term, rules) -> Term:
    child1, child2 = term1.children  
    return Term(ADD,
                [Term(MUL, [calculate_derivative(child2, rules), Term(LOG, [child1]), Term(POW, [child1, child2])]),
                 Term(MUL, [calculate_derivative(child1, rules), child2, Term(POW, [child1, Term(SUB, [child2, Term(1)])])])])


def sin_derivative(term1: Term, rules) -> Term:
    return Term(MUL, [Term(COS, [term1.children[0]]), calculate_derivative(term1.children[0], rules)])


def cos_derivative(term1: Term, rules) -> Term:
    return Term(MUL, [Term(-1), Term(SIN, [term1.children[0]]), calculate_derivative(term1.children[0], rules)])


def tan_derivative(term1: Term, rules) -> Term:
    arg = Term.copy(term1.children[0])
    return Term(ADD, [calculate_derivative(arg, rules), Term(MUL, [calculate_derivative(arg, rules), Term(POW, [term1, Term(2)])])])


def log_derivative(term1: Term, rules) -> Term:
    return Term(DIV, [calculate_derivative(term1.children[0], rules), term1.children[0]])


rules_for_differentiation = {MUL: mul_derivative, ADD: sum_derivative, SUB: sum_derivative,
                             COS: cos_derivative, SIN: sin_derivative, TAN: tan_derivative,
                             POW: power_derivative, DIV: div_derivative, LOG: log_derivative}


def derivative(term1):
    return calculate_derivative(term1, rules_for_differentiation)

MATLAB符号微分

求函数 $y=f(x{)}^2\frac{df}{dx}$ 关于 $f(x)$ 的导数

syms f(x) y
y = f(x)^2*diff(f(x),x);
Dy = diff(y,f(x))

$\begin{array}{rl}& Dy=\\ & 2f(x)\frac{\partial }{\partial x}f(x)\end{array}$

找到描述质量弹簧系统运动的欧拉-拉格朗日方程。定义系统的动能和势能。

syms x(t) m k
T = m/2*diff(x(t),t)^2;
V = k/2*x(t)^2;

定义拉格朗日量。

L = T - V

$\begin{array}{rl}& L=\\ & \frac{m{\left(\frac{\partial }{\partial t}x(t)\right)}^2}{2}-\frac{kx(t{)}^2}{2}\end{array}$

欧拉-拉格朗日方程如下

$0=\frac{d}{dt}\frac{\partial L(t,x,\dot{x})}{\partial \dot{x}}-\frac{\partial L(t,x,\dot{x})}{\partial x}$

计算项 $\partial L/\partial \dot{x}$

D1 = diff(L,diff(x(t),t))

$\begin{array}{rl}& D1=\\ & m\frac{\partial }{\partial t}x(t)\end{array}$

计算第二项 $\partial L/\partial x$

D2 = diff(L,x)

$D2(t)=-kx(t)$

找到质量弹簧系统的欧拉-拉格朗日运动方程。

diff(D1,t) - D2 == 0

$\begin{array}{rl}& \mathrm{ans}(t)=\\ & m\frac{{\partial }^2}{\partial t^2}x(t)+kx(t)=0\end{array}$

👉更新：亚图跨际