某厂每次生产时需要付生产准备费5000元，产量大于需求时多出的部分每天每件需要支付储存费1元，产品日需求量为100件，设每次产量可以任意取值，求该厂多少天生产一次，每次产量多少，可以使得总费用最小？

1.无缺货情形

首先该模型不允许缺货，因此每次产量不能低于生产周期内的需求量。其次过剩生产也没有好处，因为一个周期内的需求量是固定的，一个周期内的过剩生产导致的库存积压只会导致不必要的储存费增加。所以，每次产量必定等于生产周期内的需求量，也就是天数*100。

生产周期短，则生产所需的生产准备费高，而储存费低，生产周期长则相反。因此这是一个优化问题，我们需要建立目标函数。

评估总费用需要所有方案都有一段相同的作用时间，用天作为这个作用时间，我们只需要计算每种方案在周期内的日均费用即可。

设该厂每次生产时需要付生产准备费N元，产量大于需求时多出的部分每天每件需要支付储存费n元，产品日需求量为m件，设该厂x天生产一次，每次生产x·m件，计算一个周期内的日均费用，为： $\frac{N+(x-1)·m·n+(x-2)·m·n+...+m·n}{x}=\frac{N+\sum_{i=1}^{x-1}i·m·n}{x}=\frac{N}{x}+\frac{m·n}{x}·\frac{(x-1)x}{2}=\frac{N}{x}+\frac{mnx}{2}-\frac{mn}{2}$，对该式求导即可计算最小值。

输入N、n、m，输出导函数和使导函数为0的x，python代码实现如下：

import sympy as sp

N=float(input())
n=float(input())
m=int(input())
x = sp.Symbol('x')
fx=N/x+m*n*x/2-m*n/2
dfx=sp.diff(fx,x)
print(dfx)
x=sp.solve(dfx,x)
print(x)

输入5000、1、100时，结果如图：

绘制原函数与导函数图像：

发现在x=10时，导函数为0，原函数取最小值。

2.允许缺货情形

如果我们在一个周期内缺货，则储存费用会降低，而由于生产件数是任意的，我们只需在下次生产时加大产量以补足缺货即可。现再假设出现缺货时，每天每件缺货损失费为c，则设该厂每次生产时需要付生产准备费N元，产量大于需求时多出的部分每天每件需要支付储存费n元，产品日需求量为m件，设该厂x天生产一次，每次生产q件，则日均费用为： $\frac{N+n·\sum_{i=1}^\frac{q}{m}(q-mi)+c·\sum_{k=\frac{q}{m}}^xm(k-\frac{q}{m})}{x}$。该模型可能出现小数，为了简化，我们将求和运算变为积分运算，即 $\frac{N+n·\int_0^\frac{q}{m}(q-mi){\rm d}i+c·\int_{\frac{q}{m}}^xm(k-\frac{q}{m}){\rm d}k}{x}$。这是关于q和x的二元函数，分别对q和x求偏导，求使得偏导为0的q和x，再求二阶偏导，利用二元函数的极值判别法来计算极值，并比较获得最值。

输入N、n、m、c，输出最低日均费用和生产周期与产量的取值，python代码实现如下：

import numpy as np
import pandas as pd
import sympy
from matplotlib import pyplot

N=float(input())
n=float(input())
m=int(input())
c=float(input())

x,q,i,k= sympy.symbols('x q i k')
f=(N+n*sympy.integrate(q-m*i,(i,0,q/m))+c*sympy.integrate(m*(k-q/m),(k,q/m,x)))/x
dfx=sympy.diff(f,x)
dfq=sympy.diff(f,q)
points=sympy.solve([dfx,dfq],[x,q])
dfxx=sympy.diff(dfx,x)
dfxq=sympy.diff(dfx,q)
dfqq=sympy.diff(dfq,q)
fmin=2147483648
min_point=[]
for point in points:
    if dfxq.evalf(subs={x:point[0],q:point[1]})**2\
        -dfqq.evalf(subs={x:point[0],q:point[1]})*dfxx.evalf(subs={x:point[0],q:point[1]})<0\
             and dfxx.evalf(subs={x:point[0],q:point[1]})>0:
        if f.evalf(subs={x:point[0],q:point[1]})<fmin:
            fmin=f.evalf(subs={x:point[0],q:point[1]})
            min_point=point
print(fmin)
print(min_point)

输入为5000、1、100、2时，输出结果如图：