C#实现迭代算法：求解n级串联黑箱的最优输入x₀

问题回顾 #

给定n个串联的黑箱，每个黑箱的输入输出关系为：
$$f_i(x) = L_i \times \left( \frac{1}{\frac{x}{L_i} + \sqrt{P_i}} - \frac{1}{\sqrt{P_i}} \right)$$
输入$x_0$经过n次复合后得到$x_{n+1}=f_{n-1} \circ ... \circ f_0(x_0)$，目标是找到$x_0$最大化差值$d(x_0)=x_{n+1}-x_0$，等价于求解$d'(x_0)=0$的根。

核心推导 #

1. 单个黑箱导数：
   对$f_i(x)$求导可得：
   $$f_i'(x) = -\frac{1}{\left( \frac{x}{L_i} + \sqrt{P_i} \right)^2}$$
   二阶导数：
   $$f_i''(x) = \frac{2}{L_i \times \left( \frac{x}{L_i} + \sqrt{P_i} \right)^3}$$

2. 复合函数导数递推：
   设$\frac{dx_{k+1}}{dx_0}$为第k个黑箱输出对$x_0$的导数，递推关系为：
   $$\frac{dx_{k+1}}{dx_0} = f_k'(x_k) \times \frac{dx_k}{dx_0}$$
   初始条件$\frac{dx_0}{dx_0}=1$，最终$d'(x_0)=\frac{dx_{n+1}}{dx_0} - 1$。

迭代算法：牛顿-拉夫逊法 #

由于$d(x_0)$是凹函数（$d''(x_0)<0$），存在唯一极大值点，牛顿法可高效收敛到该点，支持任意n值：

算法步骤

1. 初始化：选择合法初始值$x_0$（如0，需满足所有黑箱输入$x_k > -L_i\sqrt{P_i}$）。
2. 正向传播：从$x_0$出发，依次计算所有黑箱的输入输出$x_1, x_2,...,x_{n+1}$，记录每个$x_k$。
3. 反向求一阶导数：从最后一个黑箱反向递推，计算$\frac{dx_{n+1}}{dx_0}$，得到$d'(x_0)$。
4. 反向求二阶导数：同样反向递推计算$d''(x_0)$，用于牛顿法的步长更新。
5. 更新$x_0$：按牛顿公式$x_0^{(k+1)}=x_0{(k)} - \frac{d'(x_0^{{(k)})}{d''(x_0}{(k)})}$更新，验证新$x_0$的合法性（避免黑箱分母无效）。
6. 收敛判断：当$|d'(x_0)| < \epsilon$或$|x_0^{(k+1)}-x_0{(k)}| < \epsilon$时停止迭代。

伪代码

FUNCTION FindOptimalX0(boxes, epsilon, maxIter):
    n = LENGTH(boxes)
    x0 = 0.0  # 初始值可根据参数调整
    FOR iter FROM 0 TO maxIter:
        # 正向计算所有x值
        x_list = [x0]
        FOR i FROM 0 TO n-1:
            L = boxes[i].L
            sqrtP = SQRT(boxes[i].P)
            denom = x_list[i] / L + sqrtP
            next_x = L * (1.0 / denom - 1.0 / sqrtP)
            x_list.APPEND(next_x)
        
        # 计算一阶导数d'(x0)
        dx = 1.0
        FOR i FROM n-1 DOWN TO 0:
            L = boxes[i].L
            sqrtP = SQRT(boxes[i].P)
            xi = x_list[i]
            denom = xi / L + sqrtP
            fi_prime = -1.0 / (denom * denom)
            dx = dx * fi_prime
        d_prime = dx - 1.0
        
        IF ABS(d_prime) < epsilon:
            RETURN x0
        
        # 计算二阶导数d''(x0)
        dx_list = [1.0]
        FOR i FROM 0 TO n-1:
            L = boxes[i].L
            sqrtP = SQRT(boxes[i].P)
            xi = x_list[i]
            denom = xi / L + sqrtP
            fi_prime = -1.0 / (denom * denom)
            dx_list.APPEND(dx_list[i] * fi_prime)
        
        current_dg = 0.0
        FOR i FROM n-1 DOWN TO 0:
            L = boxes[i].L
            sqrtP = SQRT(boxes[i].P)
            xi = x_list[i]
            denom = xi / L + sqrtP
            fi_prime = -1.0 / (denom * denom)
            fi_double_prime = 2.0 / (L * denom * denom * denom)
            current_dg = fi_double_prime * dx_list[i]^2 + fi_prime * current_dg
        d_double_prime = current_dg
        
        # 更新x0并验证合法性
        delta = d_prime / d_double_prime
        x0_new = x0 - delta
        
        # 检查新x0是否会导致黑箱输入无效
        valid = TRUE
        temp_x = x0_new
        FOR box IN boxes:
            L = box.L
            sqrtP = SQRT(box.P)
            IF temp_x <= -L * sqrtP:
                valid = FALSE
                BREAK
            temp_x = L * (1.0 / (temp_x/L + sqrtP) - 1.0/sqrtP)
        
        IF valid:
            IF ABS(x0_new - x0) < epsilon:
                RETURN x0_new
            x0 = x0_new
        ELSE:
            # 步长减半，避免无效值
            x0 = x0 - delta * 0.5
    
    RETURN x0

C#实现示例

public class BlackboxParams
{
    public double L { get; set; }
    public double P { get; set; }
}

public static double FindOptimalX0(List<BlackboxParams> boxes, double epsilon = 1e-6, int maxIter = 100)
{
    double x0 = 0.0;
    int boxCount = boxes.Count;

    for (int iter = 0; iter < maxIter; iter++)
    {
        // 正向计算所有节点的x值
        List<double> xList = new List<double> { x0 };
        foreach (var box in boxes)
        {
            double sqrtP = Math.Sqrt(box.P);
            double denom = xList.Last() / box.L + sqrtP;
            double nextX = box.L * (1.0 / denom - 1.0 / sqrtP);
            xList.Add(nextX);
        }

        // 计算一阶导数d'(x0)
        double dx = 1.0;
        for (int i = boxCount - 1; i >= 0; i--)
        {
            var box = boxes[i];
            double sqrtP = Math.Sqrt(box.P);
            double xi = xList[i];
            double denom = xi / box.L + sqrtP;
            double fiPrime = -1.0 / (denom * denom);
            dx *= fiPrime;
        }
        double dPrime = dx - 1.0;

        if (Math.Abs(dPrime) < epsilon)
            return x0;

        // 计算二阶导数d''(x0)
        List<double> dxList = new List<double> { 1.0 };
        for (int i = 0; i < boxCount; i++)
        {
            var box = boxes[i];
            double sqrtP = Math.Sqrt(box.P);
            double xi = xList[i];
            double denom = xi / box.L + sqrtP;
            double fiPrime = -1.0 / (denom * denom);
            dxList.Add(dxList[i] * fiPrime);
        }

        double currentDg = 0.0;
        for (int i = boxCount - 1; i >= 0; i--)
        {
            var box = boxes[i];
            double sqrtP = Math.Sqrt(box.P);
            double xi = xList[i];
            double denom = xi / box.L + sqrtP;
            double fiPrime = -1.0 / (denom * denom);
            double fiDoublePrime = 2.0 / (box.L * denom * denom * denom);
            currentDg = fiDoublePrime * dxList[i] * dxList[i] + fiPrime * currentDg;
        }
        double dDoublePrime = currentDg;

        // 更新x0并验证合法性
        double delta = dPrime / dDoublePrime;
        double x0New = x0 - delta;

        bool isValid = true;
        double tempX = x0New;
        foreach (var box in boxes)
        {
            double sqrtP = Math.Sqrt(box.P);
            if (tempX <= -box.L * sqrtP)
            {
                isValid = false;
                break;
            }
            tempX = box.L * (1.0 / (tempX / box.L + sqrtP) - 1.0 / sqrtP);
        }

        if (isValid)
        {
            if (Math.Abs(x0New - x0) < epsilon)
                return x0New;
            x0 = x0New;
        }
        else
        {
            x0 -= delta * 0.5;
        }
    }

    return x0;
}

关键注意点 #

• 初始值：优先选择正数，避免黑箱输入触发无效分母；若参数范围特殊，可根据$L_i$和$P_i$估算合理初始值。
• 收敛性：由于$d(x_0)$是严格凹函数，牛顿法必然收敛到唯一极大值点，无需担心局部最优问题。
• 数值稳定性：当更新后的$x_0$导致黑箱输入无效时，需调整步长（如减半），避免迭代中断。

内容的提问来源于stack exchange，提问作者Alexander Schmidt