这里有一个稍微不同的方法，它使问题保持离散，并详尽地列举了可能的解决方案。这在该情况下是可行的，因为仅存在3234个组合。在这种方法中存在对于太稀疏的数据可能已经明显的不准确性，但是显然它不排除首先拟合多项式并且然后使用那些来插值并且然后使用更细粒度的数据来抑制这种不期望的效果的可能性。在这里，每个值数组的稀疏度降低了四倍，给出了946125个组合，但运行时间仍然很短。

from itertools import product
import numpy as np

P_North = [0.11, 0.23, 0.45, 0.75, 1.11, 1.45, 1.75, 2.01, 2.17, 2.30, 2.27]
eta_North = [4.584, 18.362, 32.330, 50.288, 64.760, 74.638, 80.022, 83.494, 84.771, 85.635, 84.254]

P_West = [3.63, 4.59, 5.51, 6.50, 7.35, 8.35, 8.94]
eta_West = [82.717, 83.711, 84.341, 84.629, 83.793, 82.765, 80.756]

P_South = [4.541, 9.181, 13.771, 18.561, 22.971, 27.851, 31.821]
eta_South = [80.3471, 83.1261, 84.3591, 85.4051, 85.6741, 85.6191, 84.4601]

P_East = [18.00, 20.00, 25.00, 32.00, 35.00, 36.00]
eta_East = [81.360, 83.253, 85.589, 85.642, 85.624, 82.931]

tot_names = ["North", "West", "South", "East"]
P_tot_orig = [P_North, P_West, P_South, P_East]
eta_tot_orig = [eta_North, eta_West, eta_South, eta_East]

字符串
插值：

interp_factor = 4

def interpolate(arr, factor):
    n = len(arr) * factor
    return np.arange(min(arr) - 0.001, max(arr) + 0.001, (max(arr) - min(arr)) / n)

coefficient = [np.polyfit(P, eta, 2) for P, eta in zip(P_tot_orig, eta_tot_orig)]
polynomal = [np.poly1d(coefficient) for coefficient in coefficient]

P_tot = [list(interpolate(P, interp_factor)) for P in P_tot_orig]
eta_tot = [list(polynomal(P)) for P, polynomal in zip(P_tot, polynomal)]

型
优化：

def get_optimal_configuration(total_effect, max_relative_deviation=.05):
    min_P = total_effect
    max_P = total_effect * (1 + max_relative_deviation)
    best_E, best_P = -float('inf'), None
    powers, efficiencies = None, None
    for P_arr, eta_arr in zip(product(*P_tot), product(*eta_tot)):
        P = sum(P_arr)
        if min_P <= P <= max_P:
            E = sum(eta_arr)
            if E > best_E:
                best_P, best_E = P, E
                powers, efficiencies = P_arr, eta_arr
    if best_P is None:
        print(f'No solution for total_effect in ({min_P:.2f}, {max_P:.2f})')
    else:
        p_dist = ", ".join([f'({n}, {p:.2f})' for n, p in zip(tot_names, powers)])
        e_dist = ", ".join([f'({n}, {e:.2f})' for n, e in zip(tot_names, efficiencies)])
        print(
            f'Maximum efficiency of {best_E:.2f}% for total_effect in ({min_P:.2f}, {max_P:.2f}) reached for '
            f'total_effect = {best_P:.2f} and\npowers [{p_dist}] and\nefficiencies [{e_dist}].\n'
        )
    return best_E, best_P, powers, efficiencies

型
示例用法：

get_optimal_configuration(60)

型
印刷品

Maximum efficiency of 341.29% for P in (60.00, 63.00) reached for P = 60.58 and
powers [(North, 2.15), (West, 5.90), (South, 24.03), (East, 28.50)] and
efficiencies [(North, 84.57), (West, 84.60), (South, 85.77), (East, 86.34)].

(341.2886446785335,
 60.57911038961039,
 (2.149681818181819, 5.904714285714286, 24.025714285714283, 28.499),
 (84.57391666669983, 84.60418947837982, 85.7660321381848, 86.34450639526902))

型
说明：
首先，针对四个方向中的每一个拟合2次多项式f(P, eta)。接下来，每对数组(P, eta)通过将数组P中的点的数量粗略地增加四倍（参见interp_factor）（保持范围相同，除了一些buffer）并应用获得的多项式来获得eta的值。
然后，函数get_optimal_configuration简单地遍历每个组合，并选择任何高于指定total_effect（但不超过total_effect * max_relative_deviation）的组合，其效率之和最高（相当于最小化减去效率之和与total_effect的比率），如果存在任何这样的组合。

我不认为你做了足够的限制你的权力设置到有效的范围，在几个地方：在你的input中，以及在你的边界中。
如果您没有任何接近0的低效应值的效率数据，那么0肯定不是有效的功效效应。我首先用上界和下界进行演示。
如果强制执行这些界限，那么/4的分布显然是无效的，因为这将超出某些电厂的范围。
功率和效率的实际值的交叉比垂直线更能提供信息，并且我不认为按发电机分隔子图有多大价值-它们可以在同一个图上，特别是如果您使用半对数。

import numpy as np
import pandas as pd
from numpy.polynomial import Polynomial
from scipy.optimize import minimize, LinearConstraint
import matplotlib.pyplot as plt

df_data = (
    pd.DataFrame({
        'P': [0.11, 0.23, 0.45, 0.75, 1.11, 1.45, 1.75, 2.01, 2.17, 2.30, 2.27],
        'eta': [4.584, 18.362, 32.330, 50.288, 64.760, 74.638, 80.022, 83.494, 84.771, 85.635, 84.254],
        'loc': 'North',
    }),
    pd.DataFrame({
        'P': [3.63, 4.59, 5.51, 6.50, 7.35, 8.35, 8.94],
        'eta': [82.717, 83.711, 84.341, 84.629, 83.793, 82.765, 80.756],
        'loc': 'West',
    }),
    pd.DataFrame({
        'P': [4.541, 9.181, 13.771, 18.561, 22.971, 27.851, 31.821],
        'eta': [80.3471, 83.1261, 84.3591, 85.4051, 85.6741, 85.6191, 84.4601],
        'loc': 'South',
    }),
    pd.DataFrame({
        'P': [18.00, 20.00, 25.00, 32.00, 35.00, 36.00],
        'eta': [81.360, 83.253, 85.589, 85.642, 85.624, 82.931],
        'loc': 'East',
    }),
)
n_loc = len(df_data)
df = pd.concat(df_data)

def do_fit(group: pd.DataFrame) -> np.ndarray:
    return Polynomial.fit(x=group.P, y=group.eta, deg=2, symbol='P')

by_loc = df.groupby('loc')
min_powers = by_loc.P.min()
max_powers = by_loc.P.max()
min_effect = min_powers.sum()
max_effect = max_powers.sum()
polys = by_loc.apply(do_fit)

while True:
    try:
        total_effect = float(input(f"Power in MW this hour ({min_effect:.3f}-{max_effect:.3f}): "))
        if min_effect <= total_effect <= max_effect:
            break
    except ValueError:
        pass

def eval_polys(P: np.ndarray) -> list[float]:
    return [poly(p) for poly, p in zip(polys, P)]

def cost(P: np.ndarray) -> float:
    return -sum(eval_polys(P))

constraints = LinearConstraint(
    A=np.ones(n_loc),
    lb=total_effect, ub=total_effect,
)
bounds = np.stack((min_powers, max_powers), axis=1)

solution = minimize(fun=cost, x0=min_powers, constraints=constraints, bounds=bounds)
assert solution.success, solution.message
optimal_P = pd.Series(data=solution.x, index=polys.index)
optimal_eta = pd.Series(data=eval_polys(optimal_P), index=polys.index)
print('Optimal effect (MW):')
print(optimal_P.to_string(float_format='{:.2f}'.format))
print('Optimal eta:')
print(optimal_eta.to_string(float_format='{:.3f}'.format))

fig: plt.Figure
ax: plt.Axes
fig, ax = plt.subplots()

for loc, group in by_loc:
    P = np.linspace(start=min_powers[loc], stop=max_powers[loc], num=100)
    ax.plot(P, polys[loc](P), zorder=1, label=loc + ' (fit)')
    ax.scatter(group.P, group.eta, zorder=2, label=loc + ' (data)')

ax.scatter(optimal_P, optimal_eta, zorder=3, c='black', s=80, marker='+', label='Optimal values')

ax.set_xscale('log')
ax.set_xlabel('Effect [MW]')
ax.set_ylabel('Efficiency [%]')
ax.legend()

plt.show()

个字符
x1c 0d1x的数据
发电站“East”在P = 28.52处具有86.34%的最大效率。因此，如果我将total_effect设置为28.52，则最优解决方案将是仅运行28.52MW的East发电站，并将所有其他发电站设置为
结果可能相同，但成本标准不同：考虑到手头的实际问题，最大化效率值之和可能没有用，因为它们没有相同的权重。相反，你应该 * 最小化 * 中的总功率，根据单个效率计算。将其应用于28.52 MW的数字，并将假设更改为接受0作为有效效应，解决方案看起来更像

import numpy as np
import pandas as pd
from numpy.polynomial import Polynomial
from scipy.optimize import minimize, LinearConstraint
import matplotlib.pyplot as plt

df_data = (
    pd.DataFrame({
        'P': [0.11, 0.23, 0.45, 0.75, 1.11, 1.45, 1.75, 2.01, 2.17, 2.30, 2.27],
        'eta': [4.584, 18.362, 32.330, 50.288, 64.760, 74.638, 80.022, 83.494, 84.771, 85.635, 84.254],
        'loc': 'North',
    }),
    pd.DataFrame({
        'P': [3.63, 4.59, 5.51, 6.50, 7.35, 8.35, 8.94],
        'eta': [82.717, 83.711, 84.341, 84.629, 83.793, 82.765, 80.756],
        'loc': 'West',
    }),
    pd.DataFrame({
        'P': [4.541, 9.181, 13.771, 18.561, 22.971, 27.851, 31.821],
        'eta': [80.3471, 83.1261, 84.3591, 85.4051, 85.6741, 85.6191, 84.4601],
        'loc': 'South',
    }),
    pd.DataFrame({
        'P': [18.00, 20.00, 25.00, 32.00, 35.00, 36.00],
        'eta': [81.360, 83.253, 85.589, 85.642, 85.624, 82.931],
        'loc': 'East',
    }),
)
n_loc = len(df_data)
df = pd.concat(df_data)

def do_fit(group: pd.DataFrame) -> np.ndarray:
    return Polynomial.fit(x=group.P, y=group.eta, deg=2, symbol='P')

by_loc = df.groupby('loc')
max_powers = by_loc.P.max()
max_effect = max_powers.sum()
polys = by_loc.apply(do_fit)

while True:
    try:
        total_effect = float(input(f"Power in MW this hour (0-{max_effect:.3f}): "))
        if 0 <= total_effect <= max_effect:
            break
    except ValueError:
        pass

def eval_polys(P: np.ndarray) -> list[float]:
    return [poly(p) for poly, p in zip(polys, P)]

def power_in(Pout: np.ndarray) -> np.ndarray:
    efficiencies = eval_polys(Pout)
    return Pout / efficiencies * 100

def cost(Pout: np.ndarray) -> float:
    return power_in(Pout).sum()

constraints = LinearConstraint(
    A=np.ones(n_loc),
    lb=total_effect, ub=total_effect,
)
bounds = np.stack((np.zeros(n_loc), max_powers), axis=1)

solution = minimize(fun=cost, x0=np.zeros(n_loc), constraints=constraints, bounds=bounds)
assert solution.success, solution.message
optimal_P = pd.Series(data=solution.x, index=polys.index)
optimal_eta = pd.Series(data=eval_polys(optimal_P), index=polys.index)

print('\nOptimal effect (MW):')
print(optimal_P.to_string(float_format='{:.2f}'.format))

print('\nOptimal eta:')
print(optimal_eta.to_string(float_format='{:.3f}'.format))

fig: plt.Figure
ax: plt.Axes
fig, ax = plt.subplots()

for loc, group in by_loc:
    P = np.linspace(start=0, stop=max_powers[loc], num=100)
    ax.plot(P, polys[loc](P), zorder=1, label=loc + ' (fit)')
    ax.scatter(group.P, group.eta, zorder=2, label=loc + ' (data)')

ax.scatter(
    optimal_P.clip(lower=0.01),
    optimal_eta,
    zorder=3, c='black', s=80, marker='+', label='Optimal values')

ax.set_xscale('log')
ax.set_xlabel('Effect [MW]')
ax.set_ylabel('Efficiency [%]')
ax.legend()

plt.show()

Power in MW this hour (0-79.061): 28.52

Optimal effect (MW):
loc
East    26.46
North    2.06
South    0.00
West     0.00

Optimal eta:
loc
East    86.155
North   84.447
South   77.537
West    70.919

的字符串

的
将您建议的解决方案与上面代码中的解决方案进行比较

print(power_in(np.array([26.46, 2.06, 0, 0])).sum())
print(power_in(np.array([28.52, 0, 0, 0])).sum())

型
我们得到

33.151353481691885
33.03046970728739

型
他们很亲密。这可以通过以下方式进一步调整

• 添加雅可比矩阵（这应该足够简单），
• 使用全局优化器而不是局部优化器，
• 减少一个自由度并消除约束，

或者更好的方法是，分析计算最小值，因为（在边界内）你的成本函数是连续的和可微的。