这不是一个完整的解决方案，更像是一些可能导致一个解决方案的想法。另外，有人请再次检查我的数学。

变量和测试数据

首先，让我们从定义一些测试数据开始。在这样做的时候，我转置了a1a2a3a4矩阵，因为它会证明更方便。另外，我将它重命名为a14，因为我不想一直写整个东西，对不起。

import numpy as np
from numpy import linalg, random
from scipy import optimize

a14 = random.uniform(-1., 1., (500, 4)) * (0.25, 0.75, 1., 2.)
at = random.uniform(-1., 1., 500)
xpoints = np.cumsum(random.random(500))
initial_guess = np.full(4, 0.25)

目标
好吧，这是我第一个有点困惑的地方：OP将目标描述为两个图形之间的“面积差”。X坐标由xpoints给出。目标的Y坐标由at给出，并且我们搜索coeffs，使得a14 @ coeffs给出接近目标的Y坐标。
现在，当我们谈论“图下的面积”时，我的第一直觉是分段线性的。类似于trapezoidal rule，这给了我们公式sum{(y[i+1]+y[i]) * (x[i+1]-x[i])} / 2。我们可以稍微重新公式化，得到(y[0]*(x[1]-x[0]) + y[1]*(x[2]-x[0]) + y[2]*(x[3]-x[1]) + ... + y[n-1]*(x[n-1]-x[n-2])) / 2的形式。我们想预先计算x的差值和缩放因子。

weights = 0.5 * np.concatenate(
         ((xpoints[1] - xpoints[0],),
         xpoints[2:] - xpoints[:-2],
         (xpoints[-1] - xpoints[-2],)))

这些因子与OP使用的不同。我猜它们是对分段常数函数而不是分段线性函数进行积分？无论如何，在这两者之间切换应该足够简单，所以我将其作为替代方案留在这里。

功能丧失

根据我的解释调整OP的损失函数，得到

def loss(coeffs):
    absdiff = abs(at - (a14 * coeffs).sum(axis=1))
    return (absdiff * weights).sum()

正如我在注解中提到的，这可以通过使用矩阵向量乘法来简化。

def loss(coeffs):
    return abs(a14 @ coeffs - at) @ weights

这与abs((a14 * weights[:, np.newaxis]) @ coeffs - (at * weights)).sum()是一样的，这反过来又使我们很明显地讨论的是最小化L1范数。

a14w = a14 * weights[:, np.newaxis]
atw = at * weights
def loss(coeffs):
    return linalg.norm(a14w @ coeffs - atw, 1)

就损失函数的计算而言，这并不是一个巨大的进步，但是，它对把我们的问题变成一个规则的模式有很大的帮助。

近似解

正如在评论中所指出的，L1范数下的abs()函数是优化的毒药，因为它是不可微的。如果我们切换到L2范数，我们基本上有一个带有附加约束的最小二乘拟合。当然，这有点不合理，因为我们开始解决一个不同的，如果密切相关的问题。我怀疑解决方案将是“足够好的”;或者它可以是起始溶液，然后将其抛光以符合实际目标。
在任何情况下，我们都可以使用scipy.optimize.lsq_linear作为一个快速简单的求解器，但是它不支持线性约束，我们可以用正则化参数来模拟它。

def loss_lsqr(coeffs):
    return 0.5 * linalg.norm(a14w @ coeffs - atw)**2

grad1 = a14w.T @ a14w
grad2 = a14w.T @ at
def jac_lsqr(coeffs):
    return grad1 @ coeffs - grad2

regularization = loss_lsqr(initial_guess) # TODO: Pick better value
a14w_regular = np.append(a14w, np.full((1, 4), regularization), axis=0)
atw_regular = np.append(atw, regularization)
approx = optimize.lsq_linear(a14w_regular, atw_regular, bounds=(0., 1.))

# polish further, maybe?
bounds = ((0., 1.),) * 4
constraint = optimize.LinearConstraint(np.ones(4), 1., 1.)
second_guess = approx.x
second_guess /= linalg.norm(second_guess) # enforce constraint
exact = optimize.minimize(
        loss_lsqr, second_guess, jac=jac_lsqr,
        bounds=bounds, constraints=constraint)