你的预处理做得很差，我相信你遗漏了一个步骤。你没有展示你加载数据的方式，也没有展示你的任何框架的构造，这在这里是非常重要的。你使用了单字母变量和很多生成器，你不应该使用任何一个。作为一个起点，考虑

from functools import partial
from scipy.optimize import milp

import pandas as pd

def load() -> tuple[
    pd.DataFrame, pd.DataFrame, pd.Series, pd.DataFrame, pd.Series, pd.Series
]:
    read = partial(pd.read_excel, 'PuLP sample.xlsx')

    df_stock = read('Stock', index_col=[0, 1, 2])
    times = df_stock.columns.to_series(name='Time')
    series_stock = df_stock.stack()
    # Time is in days after Jan 1, 1900
    series_stock.index.set_names('Time', level=3, inplace=True)
    series_stock.name = 'Stock'

    series_cost = read('Cost', index_col=[0, 1]).Cost

    df_price = read('Price', index_col=[0, 1])

    # Redundant
    # series_quantity_sum = read('Custsum', index_col=[0, 1])['Sum of Quantity']

    df_order = read('Order', index_col=[0, 1, 2])

    # Redundant
    # series_cust_order = read('CustOrder', index_col=[0]).Customer

    df_cust_source = read('Customer Source')
    df_cust_source.rename(columns={'Product': 'Plant'}, inplace=True)

    return df_order, df_cust_source, series_cost, df_price, times, series_stock

def preprocess() -> pd.DataFrame:
    (
        df_order, df_cust_source, series_cost, df_price, times, series_stock,
    ) = load()

    orders_spread = pd.merge(
        left=df_order.reset_index(), right=df_cust_source, how='left',
        left_on='Customer', right_on='Customer',
    )

    orders_spread = pd.merge(
        left=orders_spread, right=series_cost, how='left',
        left_on=['Plant', 'Product'], right_on=['Plant', 'Product'],
    )

    # Some costs (per plant and product) are missing. Within one order, if a plant does not have a
    # price for all products, that plant must be eliminated.
    has_cost = (
        orders_spread
        .assign(HasCost=orders_spread.Cost.notna())
        .groupby(['Sales Order', 'Plant'])
        .HasCost.all()
    )
    orders_spread = pd.merge(
        left=orders_spread, right=has_cost[has_cost], how='inner',
        left_on=['Sales Order', 'Plant'], right_on=['Sales Order', 'Plant'],
    )

    orders_spread = pd.merge(
        left=orders_spread, right=df_price.Price, how='left',
        left_on=['Customer', 'Product'], right_on=['Customer', 'Product'],
    )

    orders_spread['Profit'] = orders_spread.Price - orders_spread.Cost

    orders_spread = pd.merge(
        left=orders_spread, right=times, how='cross',
    ).query('Time <= `Delivery Date`')

    running_stock = series_stock.groupby(level=['Plant', 'Product']).cumsum()

    orders_spread = pd.merge(
        left=orders_spread, right=running_stock,
        left_on=['Plant', 'Product', 'Time'], right_on=['Plant', 'Product', 'Time'],
    )[[
        'Sales Order', 'Product', 'Quantity', 'Plant', 'Profit', 'Time', 'Stock',
    ]].sort_values(['Sales Order', 'Product', 'Plant', 'Time'])

    return orders_spread

def main() -> None:
    # Preprocess to reduce the problem size based on order due dates and missing costs
    orders = preprocess()

    '''Maximize (price - cost), i.e. minimize (cost - price)

    1, 2, 4. Sum of quantities must equal what's asked for in the sales order
    3. Exactly one plant must satisfy each order
    5. For each day, stock taken must be <= stock avail
    6. Stock quantities in the database are cumulative
    7. Deliver on time

    There will be:
    one 'plant choice' decision variable per order per plant
        Integral and have bounds between 0 and 1.
    one 'quantity' decision variable per order per day
        Integral and have a lower bound of 0.
    
    A is a matrix of constraint rows by >28,000 decision columns.
    - For each order-product pair, there is a fixed number of ordered products.
    - For each order, there is exactly one plant.
    - For each order-day pair, the sum of stock from all plants up to that day must exceed the sum
      of all purchases up to that day
    '''

if __name__ == '__main__':
    main()

股票列只在累积求和时有用，如果使用scipy.optimize，它将变成稀疏矩阵;我不熟悉PuLP，但是有一个等价物，稍后我将详细介绍这个例子。

这里你忽略的线性规划概念是一个all or nothing约束。我没有深入到你的模型中，另一个解决方案在导入数据方面有一些好的想法。但是作为一个开始，我会做一个小得多的第二个“玩具”数据集来排除故障。
因此，您只想从一个地点采购订单，并且只有在他们有库存的情况下。因此，您需要一个二进制变量来指示订单是否从给定的工厂完成。类似于：

fulfill[order, plant] ∈ {0, 1}   # 1 if order is fulfilled from this plant

然后（如果在其他地方没有这样做），您希望通过对所有工厂的每个订单进行求和来防止多个订单履行，类似于：

for order in orders:
    sum(fulfill[order, plant] for plant in plants) <= 1

然后将所有订单的数量链接到此订单，类似于（您需要使用其他变量展开）：

for order in order:
    for plant in plants:
        deliver[order, plant] <= fulfill[order, plant] * order_quantity

然后确保您没有强制执行所有订单的约束，这可能是不可行的
然后整理你的目标函数，（1）对所有完成的订单给予信贷，（2）对成本收取一些（负的，小的）惩罚，并使其最大化。
在计算成本惩罚时，需要做一些数学运算。假设您总是在可能的情况下交付订单，您需要做的是将最大成本惩罚加权为略小于履行订单的价值，以确保求解程序总是选择交付，而不是最小化成本。