尼克提出的主要问题是正确的（尽管不全面）：有许多局部极小值，目标是不可微的。
但同样重要的是：你的目标根本不适合作为一个目标--它应该是一个约束;整个系统这意味着现在是整数线性规划的时候了，它比一般的minimize有更好的收敛性和准确性保证。（无论您选择哪种方法）。更具体地说，您正在进行约束编程，这里根本没有目标（或者你可能知道更好的目标应该是什么，比如最小化x的系数）。
除了x、w、a、“距离元素”和m的变量之外，还需要一个二进制选择矩阵，其x的倍数等于m。这通过使用所谓的大M约束来解决x乘以某个倍数得到m是非线性的问题。
我已经尝试让它易于理解，在任何时候你都可以调用some_constraint.A.toarray()并在PyCharm的矩阵可视化器中查看结果，看看发生了什么;但你可能还需要深入研究一些关于线性规划的介绍性材料，以更好地掌握它。

import numpy as np
import scipy.sparse as sp
from scipy.optimize import milp, Bounds, LinearConstraint

s = 12_500
spans = np.array((25_450, 25_100, 25_450))
ejs = np.zeros(len(spans) - 1)  # expansion or separation joints
min_bw, max_bw = 0, 1200  # light pole barrier width
consts = spans.size*2 - 1  # Number of distance elements [b,c,d...]
m_count = spans.size*2 + 1
n = (spans - s/2)//s  # consecutive spacings of barriers for a given span

# In the old solution's objective for the initial guess, m/x ranges from 3.6 to 7.3. Ballpark 20 as a maximum.
mul_count = 20
multiples = np.arange(1, 1 + mul_count)[:, np.newaxis]

'''
Variables: x, w, a, consts, m, and then
a m_count * mul_count matrix of selectors of integer multiples of m.
Each selector is binary (0 or 1), and indicates times 1, times 2, ... times mul_count respectively.
'''

# Constraint matrix for all m equations, starting with the 'm' term itself
m_A = sp.hstack((
    sp.csc_array((m_count, 3 + consts)),
    -sp.eye(m_count),
    sp.csc_array((m_count, m_count*mul_count)),
), format='lil')

# Column positions: x, w, a, consts, m, mi
# First row: a - w/2 == m0: -0.5w + a - m == 0
m_A[0, 1:3] = -0.5, 1
# Middle rows
m_A[1:, 1] = -1  # -w (middle and last row)
m_A[1:-1, 3: 3+consts] = -sp.eye(consts)  # -consts
m_A[1:-1, 3+consts: 3+consts+m_count - 2] -= sp.eye(m_count - 2)  # -m[i - 1]

m_rhs = np.empty(m_count)
# First row
m_rhs[0] = 0
# Even rows: m[i] = s - w - m[i-1] - ejs[(i - 1)//2] - consts
#            -s + ejs[(i - 1)//2] = -w - m[i] - m[i-1] - consts
m_rhs[2:-1:2] = ejs - s
# Odd rows: m[i] = spans[(i - 1)//2] - (m[i - 1] + n[(i - 1)//2] * s + w + arr[2 + i])
#            -spans[(i - 1)//2] + n[(i - 1)//2]*s = -m[i] - m[i - 1] - w - consts
m_rhs[1:-1:2] = n*s - spans
# Last row: m[-1] = s - w
#           -w - m[-1] = -s
m_rhs[-1] = -s

constraint_m = LinearConstraint(A=m_A, lb=m_rhs, ub=m_rhs)
# To see the constraint matrix: view constraint_m.A.toarray() in a matrix viewer

# Exactly one selector per m-variable must be on
constraint_selection = LinearConstraint(
    A=sp.hstack((
        sp.csc_array((m_count, 3 + consts + m_count)),
        sp.kron(
            A=sp.eye(m_count),
            B=np.ones(mul_count),
            format='csc',
        )
    )),
    lb=np.ones(m_count),
    ub=np.ones(m_count),
)

'''
Enforce the multiple constraint:
m >= x*i                     m <= x*i       (if mi is selected)
-i*x + m >= 0                i*x - m >= 0   (if mi is selected)
-i*x + m + M(1 - mi) >= 0    i*x - m + M(1 - mi) >= 0
-i*x + m - M*mi >= -M        i*x - m - M*mi >= -M
Column positions: x, w, a, consts, m, mi
'''
ix_m = sp.hstack((
    np.tile(multiples, (m_count, 1)),
    sp.csc_array((m_count * mul_count, 2 + consts)),
    sp.kron(
        -sp.eye(m_count),
        np.ones((mul_count, 1)),
    ),
))
M = 1e6  # big-M constraint coefficient
mmi = -M*sp.eye(m_count * mul_count)
constraint_mul = LinearConstraint(
    A=sp.bmat((
        ( ix_m, mmi),
        (-ix_m, mmi),
    )),
    lb=np.full(shape=2*m_count*mul_count, fill_value=-M),
    ub=np.inf,
)

result = milp(
    c=np.zeros(3 + consts + m_count + m_count*mul_count),  # There is no optimization objective
    integrality=np.concatenate((
        np.zeros(3 + consts + m_count),  # x, w, a, consts and m continuous
        np.ones(m_count * mul_count),    # mi integral
    )),
    bounds=Bounds(
        lb=np.concatenate((
            (1_500, 1_500, 3_000),  # x, w, a
            np.full(shape=consts, fill_value=min_bw),  # consts
            np.zeros(m_count + m_count*mul_count),     # m, mi
        )),
        ub=np.concatenate((
            (2_100, 2_600, s),  # x, w, a
            np.full(shape=consts, fill_value=max_bw),   # consts
            np.full(shape=m_count, fill_value=np.inf),  # m unbounded
            np.ones(m_count*mul_count),  # mi are binary, so 0-1
        )),
    ),
    constraints=(
        constraint_m,
        constraint_selection,
        constraint_mul,
    ),
)
assert result.success, result.message
(x, w, a), distances, m, mi = np.split(result.x, (3, 3+consts, 3+consts+m_count))
mi = mi.reshape((m_count, mul_count))
multipliers = (mi @ multiples).ravel()

print(result.message)
print(f'x = {x}')
print(f'w = {w}')
print(f'a = {a}')
print(f'distance elements = {distances}')
print(f'm = {m}')
print(f'm/x quotients: {m/x}')
print(f'multipliers = {multipliers}')
print('multiply selectors =')
print(mi)

个字符
这个问题的规模相当小，而且稀疏性很高，因此可以快速求解。

plt.imshow(
    np.sign(
        sp.vstack((
            constraint_m.A,
            constraint_selection.A,
            constraint_mul.A,
        )).toarray()
    ),
)
plt.axis('scaled')
plt.show()

型
x1c 0d1x的数据

你现在有两个问题。

1. m % x有一个使优化器混淆的导数。
   1.你要解决的问题有很多局部最小值。
   让我们首先为优化器提供有关如何优化函数的更多信息。

混淆导数

函数minimize()是一个局部优化器，如果有另一个最小值，它必须爬上一座山才能到达，它就不能到达它。
假设你想最小化函数f（x）= x % 100，其中x限制在50和150之间。
x1c 0d1x的数据
假设优化器从x = 90开始，它查看导数，然后沿着斜率向左移动，然后在x = 50时卡住，同样的事情也发生在你的问题中。
要解决这个问题，你可以创建一个函数，奖励优化器接近x = 100。我想出了这个函数，它是一个余弦波，当variable接近target的倍数时，它接近于零。

def modular_difference(variable, target):
    # Move valley very slightly forward
    overshoot_distances_amount = 1e-8
    distance_to_next = variable % target / target
    distance_to_next -= overshoot_distances_amount
    ret = -np.cos(2*np.pi*distance_to_next) + 1
    return ret

字符串
(Note：overshoot_distances_amount的用途将在后面解释。）
这给出了下图：

的
这告诉优化器，如果它在x=99，它应该向x=100移动。
我用下面的方法将这个新函数集成到您的代码中。
首先，我改变了原始的目标函数，使其返回m和x，称之为objective_inner()。

def objective_inner(arr, n, s, spans, ejs):
    x = arr[0]
    w = arr[1]
    a =  arr[2] # (the offset is until the mid axis of the light pole)
    m = np.ones(len(spans) * 2 + 1) # We make room for the last constraint h=(s-w)
    m[0] = a - w / 2 

    for i in range(1, len(spans) * 2):
        j = (i - 1) // 2
        if i % 2 == 0:
            m[i] = (s - w) - (m[i - 1] + ejs[j] + arr[2 + i])
        else:
            m[i] = spans[j] - (m[i - 1] + n[j] * s + w + arr[2 + i])
    m[-1] = s - w  # We also need the light distance between two poles to be divisible by x
    return m, x

型
然后，我写了两个版本的目标函数。

def objective_differentiable(arr, n, s, spans, ejs):
    m, x = objective_inner(arr, n, s, spans, ejs)
    return modular_difference(m, x).sum()

def objective_old(arr, n, s, spans, ejs):
    m, x = objective_inner(arr, n, s, spans, ejs)
    return np.sum(m%x)

型
第一个函数是你的目标的一个版本，它有很好的导数，使用正弦波。第二个函数和你原来的目标一样。
然而，修正导数本身是不够的，它仍然会陷入局部极小值。

全局优化器

有些函数不能被局部优化器优化，因为优化器会陷入局部最小值。这就是其中之一。
为了解决这个问题，我使用了basinhopping函数，它使用了一个局部优化器，但也在随机方向上采取了很大的步骤，希望找到一个更好的最小值。
我发现下面的工作很好。

minimizer_kwargs = dict(bounds=all_bounds, method='Nelder-Mead', args=(n, s, spans, ejs))
sol = basinhopping(objective_differentiable, x0=initial, minimizer_kwargs=minimizer_kwargs, niter=100)
print("basinhopping sol", sol)
sol = minimize(fun=objective_old, x0=sol.x, bounds=all_bounds, method='Nelder-Mead', args=(n, s, spans, ejs))
print("minimize true sol", sol)

型
这将执行以下操作：

• 首先，它使用带有Nelder-Mead的basinhopping作为局部优化器，优化目标的可微版本。
• 然后，它将找到的解传递给你的原始目标。我这样做是为了比较它是否找到了比原始目标更好的解。（注意：这就是为什么上面的代码中存在overshoot_distances_amount。这会选择一个稍微在模循环点之后的解，这使得优化器更容易找到objective_old()的精确最小值。）

这可以找到在原始目标上获得1 e-4或更低分数的解决方案。
以下是它找到的解决方案之一：

array([1501.43791458, 1989.93457535, 5499.28104559,  449.99999353,
          0.000001  ,  100.00001273,    0.00000361,  450.00000564])

型