这是困难的，我不相信我的解决方案是最有效的;但它基本上是有效的。大致来说
1.在径向表面pi的中心内，找到其最小值的位置。
1.开始一个广泛的搜索，引用了最小的位置，任何对低错误，使用一个调用，以最小化。
1.从该对的原点开始，使用有界跟随例程向左和向右进行分支。如果以下任何一项变为真，则终止：

• 曲线出界。
• 组件之间的误差超过公差。
• 对收敛到一个临界点。

from typing import NamedTuple

import pandas as pd
import scipy.optimize
import numpy as np
import matplotlib.pyplot as plt

class Params(NamedTuple):
    uaa: float
    ubb: float
    uab: float
    na: int
    nb: int

def F(
    phi: np.ndarray,
    uaa: float, ubb: float, uab: float, na: int, nb: int,
) -> np.ndarray:
    """
    Helmholtz-free energy for a Flory-Huggins-like model

    :param phi: (2xAxB) array of phi independent coordinates
    :return: AxB array of function values
    """
    coef_shape = (2,) + (1,)*(len(phi.shape) - 1)
    nanb = np.array((na, nb)).reshape(coef_shape)
    uaub = np.array((uaa, ubb)).reshape(coef_shape)
    tot = phi.sum(axis=0)

    entropy = (
        (phi/nanb * np.log(phi)).sum(axis=0)
        + (1 - tot)*np.log(1 - tot)
    )
    energy = (
        phi**2 * (-0.5 * uaub)
    ).sum(axis=0) - uab * phi.prod(axis=0)
    return entropy + energy

def gradF(
    phi: np.ndarray,
    uaa: float, ubb: float, uab: float, na: int, nb: int,
) -> np.ndarray:
    """
    First-order Jacobian of F in phi.
    :param phi: (2xAxB) array of phi independent coordinates
    :return: (2xAxB) array of gradient vectors
    """
    coef_shape = (2,) + (1,)*(len(phi.shape) - 1)
    nanb = np.array((na, nb)).reshape(coef_shape)
    uaub = np.array((uaa, ubb)).reshape(coef_shape)
    return (
        (np.log(phi) + 1) / nanb
        - np.log(1 - phi.sum(axis=0))
        - phi * uaub
        - phi[::-1] * uab
        - 1
    )

def hessianF(
    phi: np.ndarray,
    uaa: float, ubb: float, uab: float, na: int, nb: int,
) -> np.ndarray:
    """
    Hessian of F in phi.
    :param phi: (2xAxB) array of phi independent coordinates
    :return: (2x2xAxB) array of Hessian matrices
    """
    coef_shape = (2,) + (1,)*(len(phi.shape) - 1)
    nanb = np.array((na, nb)).reshape(coef_shape)
    uaub = np.array((uaa, ubb)).reshape(coef_shape)

    diag = 1/phi/nanb - uaub
    antidiag = -uab
    hessian = np.empty(shape=(2, 2, *phi.shape[1:]), dtype=phi.dtype)
    hessian[[0, 1], [0, 1], ...] = diag
    hessian[[0, 1], [1, 0], ...] = antidiag
    hessian += 1/(1 - phi.sum(axis=0))

    return hessian

def test_grad(params: Params) -> None:
    """
    Sanity test to ensure that analytic gradients are correct.
    """
    xk = np.array((0.2, 0.3))
    assert scipy.optimize.check_grad(F, gradF, xk, *params) < 1e-6

    estimated = scipy.optimize.approx_fprime(xk, gradF, 1e-9, *params)
    analytic = hessianF(xk, *params)
    assert np.allclose(estimated, analytic, rtol=0, atol=1e-6)

def mu(
    phi: np.ndarray,
    uaa: float, ubb: float, uab: float, na: int, nb: int,
) -> np.ndarray:
    """
    Chemical potential (first-order Jacobian of F)
    """
    return gradF(phi, uaa, ubb, uab, na, nb)

def Pi(phi: np.ndarray, params: Params) -> np.ndarray:
    """Osmotic pressure"""
    mu_ab = mu(phi, *params)
    return (mu_ab * phi).sum(axis=0) - F(phi, *params)

def plot_isodims(
    phi: np.ndarray,
    mua: np.ndarray,
    mub: np.ndarray,
    pi: np.ndarray,
) -> plt.Figure:
    """
    Contour plot of individual mua/b/pi series, all over phi_a and phi_b
    """
    fig: plt.Figure
    axes: tuple[plt.Axes, ...]
    fig, axes = plt.subplots(nrows=1, ncols=3)
    fig.suptitle('Individual series')

    (ax_mua, ax_mub, ax_pi) = axes
    for ax in axes:
        ax.set_xlabel('phi a')
        ax.grid(True)
    ax_mua.set_ylabel('phi b')
    ax_mub.set_yticklabels(())
    ax_pi.set_yticklabels(())

    ax_mua.contour(*phi, mua)
    ax_mua.set_title('mu a')
    ax_mub.contour(*phi, mub)
    ax_mub.set_title('mu b')
    ax_pi.contour(*phi, pi)
    ax_pi.set_title('pi')

    return fig

def plot_example_proximity(
    phi: np.ndarray,  # 2x500x500
    mua: np.ndarray,  # 500x500
    mub: np.ndarray,  # 500x500
    pi: np.ndarray,   # 500x500
    params: Params,
    phi_pimin: np.ndarray,      # 2x
    phi_endpoints: np.ndarray,  # 2x2
    zapprox: np.ndarray,        # 500x500
    followed_phi: np.ndarray,   # 2x2x343-ish
) -> plt.Figure:
    """
    Plot, for the initial-fit origin point (1) and its pair (2), all three isocurves
    Include the Hessian determinant estimation (zapprox)
    """
    phi0 = phi_endpoints[:, 0]  # origin point (1)
    phi1 = phi_endpoints[:, 1]  # next point (2)
    mua0, mub0 = mu(phi0, *params)   # mu values at origin
    pi0 = Pi(phi0, params)  # pi value at origin

    # Array indices where mu closely matches that at the origin
    mua_close_y, mua_close_x = (np.abs(mua - mua0) < 1e-3).nonzero()
    mub_close_y, mub_close_x = (np.abs(mub - mub0) < 1e-3).nonzero()

    # Error from the pi value seen at the origin
    pi_close = pi - pi0

    fig: plt.Figure
    ax: plt.Axes
    fig, ax = plt.subplots()
    ax.set_title('Isocurves from initial fit origin to next,'
                 '\nwith Hessian determinant')

    ax.scatter(*phi0, c='black', marker='+', s=100, label='origin (1)')
    ax.scatter(*phi1, c='pink', marker='+', s=100, label='next (2)')
    ax.scatter(*phi_pimin, c='blue', marker='+', s=100, label='pi_min')
    ax.plot(phi[0, 0, mua_close_x], phi[0, 0, mua_close_y], label='mua')
    ax.plot(phi[1, mub_close_x, 0], phi[1, mub_close_y, 0], label='mub')
    ax.contour(*phi, pi_close, levels=[0], colors='purple')  # Pi isocurve
    ax.contour(*phi, zapprox, levels=[0], colors='brown')  # Hessian determinant
    ax.plot(*followed_phi[:, 0, :], label='follow1')
    ax.plot(*followed_phi[:, 1, :], label='follow2')

    ax.legend()

    return fig

def pairwise_error(phi: np.ndarray, params: Params) -> float:
    """
    For a phi pair, calculate the error between pair components.
    :param phi: 2x2, a1 a2 b1 b2
    :return: Least-squares error between mu_a, mu_b and pi
    """
    mua, mub = mu(phi, *params)
    pi = Pi(phi, params)
    return (
        (mua[0] - mua[1])**2 +
        (mub[0] - mub[1])**2 +
        10 * (pi[0] - pi[1])**2
    )

def least_sq_difference(phi: np.ndarray, params: Params) -> float:
    """
    From given phi endpoints, calculate exact (not estimated) values for mu and pi, and return
    their least-squared distance as the objective cost
    """
    phi = phi.reshape((2, 2))
    return pairwise_error(phi, params)

def initial_pair_estimate(
    phi_min: float, phi_max: float,
    phi_pimin: np.ndarray,  # 2x array
    params: Params,
) -> np.ndarray:
    """
    Perform an initial fit to find any pair on the target curve

    :param phi_min: Minimum search space bound of phi in both axes
    :param phi_max: Maximum search space bound of phi in both axes
    :param phi_pimin: Phi coordinates of minimum of pi; the search space is centred here
    :param params: uaa, ubb, uab, na, nb
    :return: endpoint coordinates in phi (2x2)
    """
    phia_pimin, phib_pimin = phi_pimin

    # As in original problem construction, the first point must be below-left of the second point;
    # and the distance must be greater than 0 to avoid degeneracy (superimposed pair).
    nondegenerate = scipy.optimize.LinearConstraint(
        A=[
            [-1, 1,  0, 0],
            [ 0, 0, -1, 1],
        ],
        lb=0.05,
    )

    result = scipy.optimize.minimize(
        fun=least_sq_difference, args=(params,),
        x0=(
            phia_pimin,   phia_pimin*3/2,
            phib_pimin/2, phib_pimin,
        ),
        # The origin is below-left of the minimum location of pi
        bounds=scipy.optimize.Bounds(
            lb=[phi_min,        phi_min,    phi_min,    phi_min],
            ub=[phia_pimin*0.5, phi_max,    phib_pimin, phi_max],
        ),
        constraints=nondegenerate,
        tol=1e-12,
    )
    if not result.success:
        raise ValueError(result.message)

    endpoints = result.x.reshape((2, 2))
    print(f'Initial pair found with error={result.fun:.2e}:\n{endpoints}')
    return endpoints

def characterise_initial(
    params: Params,
    phi_min: float = 0,
    phi_max: float = 1,
) -> tuple[
    np.ndarray,  # mesh-like phi coordinate space, 2xAxB
    np.ndarray,  # mua, AxB
    np.ndarray,  # mub, AxB
    np.ndarray,  # pi, AxB
    np.ndarray,  # Hessian determinant estimate, AxB
]:
    """
    Fill out a bunch of initial dimension information
    """
    phi_a = phi_b = np.linspace(phi_min, phi_max, 500)
    phi = np.stack(np.meshgrid(phi_a, phi_b))
    mua, mub = mu(phi, *params)
    pi = Pi(phi, params)

    hess = hessianF(phi, *params)
    zapprox = np.linalg.det(hess.transpose(2, 3, 0, 1))

    return phi, mua, mub, pi, zapprox

def estimate_pimin(phi: np.ndarray, pi: np.ndarray) -> np.ndarray:
    """
    Start referenced from the (estimated) minimum of pi, in the middle of the region of interest
    The Hessian estimate runs through this point.
    """
    ijmin = np.unravel_index(pi.argmin(), pi.shape)
    phi_pimin = phi[:, ijmin[0], ijmin[1]]
    print(f'Pi minimum point of {pi[ijmin]:.5f} at {phi_pimin}')
    return phi_pimin

def update_bounds(
    phi_old: np.ndarray,
    phi_new: np.ndarray,
    step: float,
) -> tuple[
    np.ndarray,  # new solution estimate x0
    np.ndarray,  # new lower bounds
    np.ndarray,  # new upper bounds
]:
    """
    After a follow-solve step, move the bounds. These bounds constitute a small moving rectangle
    with some wiggle room centered on the next place we guess a solution point will appear (using
    linear extrapolation). This makes it so the inner minimizer doesn't need to work as hard.
    :param phi_old: Previous phi, 2x2
    :param phi_new: New phi, 2x2
    :param step: Approximate inter-point distance
    :return: Updated guess for next solution, and bounds for next minimize call
    """
    delta = phi_new - phi_old

    # Determine which is the major axis of motion; assign the other to have more freedom
    if abs(delta[0, 0]) > abs(delta[0, 1]):
        wide_axis = 1
    else:
        wide_axis = 0

    # Put the next guess roughly at distance 'step' in direction 'delta'
    scale = step / np.sqrt(delta[0, 0]**2 + delta[1, 0]**2)
    x0 = phi_new + delta * scale

    # One new bound is halfway between points 'new' and 'x0'. The other is past x0.
    boundp0 = (x0[:, 0] + phi_new[:, 0])/2
    boundq0 = boundp0 + delta[:, 0]*scale

    lobound = np.array((
        (min(boundp0[0], boundq0[0]), x0[0, 1] - 10*step),
        (min(boundp0[1], boundq0[1]), x0[1, 1] - 10*step),
    ))
    hibound = np.array((
        (max(boundp0[0], boundq0[0]), x0[0, 1] + 10*step),
        (max(boundp0[1], boundq0[1]), x0[1, 1] + 10*step),
    ))

    # For the off-axis (not the major axis of motion), bound width is wider
    lobound[wide_axis, 0] -= 10*step
    hibound[wide_axis, 0] += 10*step

    # Keep step bounds within overall bounds
    lobound = lobound.clip(min=0)
    hibound = hibound.clip(max=1)

    return x0, lobound, hibound

def check_critical(
    direction: int,
    phi: np.ndarray,
    step: float,
    params: Params,
) -> bool:
    """
    If the two solution points have converged on one location, this is a critical point; report it
    and tell the outer optimizer to terminate this branch.
    :param direction: 1 or -1 (right or left)
    :param phi: New solution point, 2x2
    :param step: Approximate inter-point distance
    :param params: Passed to the chemical routines
    :return: True if this is a critical point.
    """
    if not np.all(
        np.abs(phi[:, 0] - phi[:, 1]) < step
    ):
        return False

    print(f'Convergence in direction {direction}: critical point at')
    df = pd.DataFrame(
        {
            'phia': phi[0],
            'phib': phi[1],
            'mua': mu(phi, *params)[0],
            'mub': mu(phi, *params)[1],
            'pi': Pi(phi, params),
        },
    )
    print(df.T)
    return True

def follow(
    phi_endpoints: np.ndarray,
    params: Params,
    step: float = 1e-3,
    tol: float = 1e-9,
    maxiter: int = 1000,
) -> np.ndarray:
    """
    Multidimensional, stepped following algorithm.
    :param phi_endpoints: Initial phi search point; will branch in two directions from here
    :param params: Passed to chemical routines
    :param step: Roughly, distance between points. Actual distance will vary
    :param tol: If we follow the curve to a point where the error between components exceeds this
                tolerance, we bail.
    :param maxiter: Upper limit on the number of points per branch
    :return: Array of phi values, separated by branch.
    """

    phi_halves = []
    error_min = np.inf
    error_max = -np.inf

    # Search left and right from initial-fit origin
    for direction in (-1, +1):
        phi_old = phi_endpoints.copy()
        phi_built = [phi_old]
        phi_halves.append(phi_built)
        off0 = direction*step
        off1 = off0/2
        lobound = phi_endpoints + (
            (min(off0, off1), -step),
            (-step,           -step),
        )
        hibound = phi_endpoints + (
            (max(off0, off1), step),
            (step,            step),
        )
        x0 = (lobound + hibound)/2

        for _ in range(maxiter):
            result = scipy.optimize.minimize(
                fun=least_sq_difference, args=(params,),
                x0=x0.ravel(), tol=1e-12,
                bounds=scipy.optimize.Bounds(lb=lobound.ravel(), ub=hibound.ravel()),
            )
            phi_new = result.x.reshape(2,2)
            if not result.success:
                # This is not fatal in many cases.
                if result.message != 'ABNORMAL_TERMINATION_IN_LNSRCH':
                    print(result.message)
                    break
            if result.fun > tol:
                print(f'Search direction {direction}: out of tol; terminating')
                break

            error_max = max(result.fun, error_max)
            error_min = min(result.fun, error_min)
            phi_built.append(phi_new)

            if check_critical(direction, phi_new, step, params):
                break

            x0, lobound, hibound = update_bounds(phi_old, phi_new, step)
            if np.any(hibound - lobound <= 0):
                print(f'Search direction {direction}: out of bounds; terminating')
                break

            phi_old = phi_new

    print(f'Follow complete; error within {error_min:.2e} and {error_max:.2e}')

    phi_series = np.stack(
        phi_halves[0][::-1] + phi_halves[1],
        axis=-1,
    )
    return phi_series

def main() -> None:
    params = Params(uaa=0, ubb=0, uab=4.36, na=10, nb=6)
    phi_min, phi_max = 1e-6, 0.4

    test_grad(params)

    phi, mua, mub, pi, zapprox = characterise_initial(params, phi_min, phi_max)
    phi_pimin = estimate_pimin(phi, pi)
    phi_endpoints = initial_pair_estimate(phi_min, phi_max, phi_pimin, params)

    followed_phi = follow(phi_endpoints, params)

    # plot_isodims(phi, mua, mub, pi)
    plot_example_proximity(
        phi, mua, mub, pi, params, phi_pimin,
        phi_endpoints, zapprox, followed_phi,
    )
    plt.show()

if __name__ == '__main__':
    main()

Pi minimum point of -0.01350 at [0.24448937 0.22925894]
Initial pair found with error=7.35e-13:
[[0.12118504 0.35659973]
 [0.0191144  0.19124418]]
Convergence in direction -1: critical point at
             0         1
phia  0.058309  0.058303
phib  0.340558  0.340551
mua  -2.160096 -2.160096
mub  -0.758148 -0.758148
pi    0.086084  0.086084
Convergence in direction 1: critical point at
             0         1
phia  0.302626  0.302102
phib  0.069781  0.069417
mua  -0.857908 -0.857908
mub  -2.130651 -2.130651
pi    0.043277  0.043276
Follow complete; error within 3.12e-17 and 4.33e-10

所有三个组件的简单图来观察它们是如何工作的：

这个比较复杂：

• 黑十字：初始拟合，原点对的一半
• 粉红色的十字架：初始配合，下半对
• 蓝十字：最小pi的位置
• 橙子和蓝色线：初始拟合对mu_a和mu_B的等曲线
• 紫色曲线：初始拟合对的等曲线，pi
• 布朗环：Hessian行列式
• 绿色：起点-
• 红色：曲线的后半部分，从交叉分支向左和向右延伸到与绿色汇合的临界点

在最近的编辑中，修复了剩余的bug：

• 用于跟随解决方案的步进边界需要是长条形的，在与运动主轴正交的轴上具有宽公差
• 海森近似有一个换位错误，使事物显得非物理