制砖机生产计划.docx

制砖机生产计划摘要某重型机械厂给出了其2011年上半年的销售预测情况，要求在成本波动的情况下制定出上半年的生产计划，使利润最大。本文针对是否对该厂月生产能力进行限制，分别构建相应的目标函数和约束条件，并通过转化为动态规划模型结合MATLAB编程，求得不同条件下的最优生产计划。问题一：在生产能力没有限制的情况下，考虑到销售额、固定成本和销售费用不随生产数量和建材价格的波动而波动，本文将目标函数简化为最小化变动成本与贮存费之和，并构建相应约束条件，然后转化为动态规划模型，结合MATLAB编程进行求解，简化了求解过程。最终该厂每月的生产台数分别为：33，34，31，31，29，29。问题二：在问题一的基础上，对每个月的生产能力进行限制。本文在模型一的基础上，增加了一个约束条件，并相应的更改了MATLAB中的M函数文件，得到添加约束条件后的最优生产计划。最终该厂每月的生产台数分别为：33，33，31，32，29，29。关键词：最小成本多阶段决策动态规划MATLAB1、问题提出的背景2010年以来，央行和发改委出台了一系列措施平抑建材价格，但由于对建材需求结构而言，总体上求大于供的市场状况没有得到根本改善，预计2011年的建材价格仍会有一定的增长。因此，在需要建材的制造业中，需要按照建材价格的增长幅度结合贮存成本等其他费用合理地控制生产计划，以达到利润最大化。现某重型机械厂通过对历史资料进行回归分析（即数据拟合），并给合今年上半年可能出现的影响制砖机销售的因素，预测该厂2011年上半年的销售情况如表1.1所示：月份123456增长速度10%10%20%20%30%30%月份123456销售量（台）423241672529表1.12011上半年销售预测表该厂的制砖机2010年12月的销售均价为48万元/台，今年上半年的售价保持不变。2010年12月末尚有49台未售出。制砖机从计划生产到售出会发生下列费用：（1）生产成本，包括固定成本（主要是指厂房、机器设备的折旧）和可变成本（钢材、其他材料和人工成本等，其中人工成本在可变成本中占到大约40%），按照2010年12月份的建材价格计算，可变成本（万元）与制砖机生产台数的平方成正比，比例系数为0.5。且可变成本与建材价格上涨幅度有关，例如建材价格上涨10%，则可变成本是按前面方法计算结果的1.1倍。（2）销售费用，与当月的销售金额成正比。（3）贮存费，生产出的制砖机未售出的必须贮存，即该厂生产的制砖机平均每台每月的贮存费为0.1万元。预计今年上半年建材的价格仍会有一定的增长。预计的增长速度（以2010年12月的价格为基准）见表1.2：表1.22011上半年建材价格增长表2、问题的重述与分析2.1问题的重述该重型机械厂打算在上半年将生产的制砖机全部售完，希望制定一个合理的生产计划在成本波动的情况下得到最大利润。其中，成本包括生产成本、销售费用和贮存费，生产成本随建材价格的波动而起伏。问题一：如果该厂的月生产能力没有限制，并且允许期货销售，如何制定月生产计划？问题二：如果该厂每月的生产能力限于33台，并且允许期货销售，该如何制定月生产计划？2.2问题的分析（1）对问题一的分析为得到该厂最优生产计划，需建立一个以最大利润为目标的模型，由于：=-利润单价销售量生产成本销售费用贮存费（2-1）其中：生产成本为固定成本和可变成本。固定成本为定值，不受制砖机生产台数和价格波动的影响，可变成本为关于制砖机生产台数和建材价格上涨幅度的函数；销售费用与当月销售额成正比，故也为定值，不受制砖机生产台数和价格波动的影响；贮存费与制砖机库存有关。由于固定成本、销售额和销售费用都为定值，不受制砖机生产台数和价格波动的影响，所以最终建立的目标函数为：min+(可变成本贮存费)（2-2）列出的约束条件如下：月初库存约束。本月初库存为上月初库存量与上月制造量之和，再减去上月销售量得到的差值。月初库存允许为负值，此时的数值代表上月的期货交易量；初始值。第一个月的库存为2010年12月的库存量，题目中已知，为49台；非负约束。每月生产的量不为负数。在求解该模型时，本文首先将模型转化为动态规划模型，然后用MATLAB编程进行求解。（2）对问题二的分析问题二在问题一的基础上，对每个月的生产能力进行限制，本文在模型一基础上增加了一个生产量的限制条件，相应的对动态模型进行改进，最终对MATLAB中的M函数文件进行更改，计算得到最终该厂的生产计划。3、问题的假设（1）假设库存足够大，不考虑库存容量问题；（2）假设月销售量预测准确，不考虑因实际月销售量波动导致的制砖机生产计划波动；（3）假设不同时段生产的制砖机的实际价值一样，不考虑库存制砖机的折旧率；（4）假设月初开始生产，月末进行销售；（5）假设该工厂生产能力足够大；4、符号说明i：上半年月份的编号（其中i=1,2,3,4,5,6）；p0：制砖机的销售价格（根据题目所给的信息可知p0是个常数）；xi：第i个月的生产计划；a0：该厂上半年的固定成本；ai：第i个月建材价格的增长速度；k：销售费用与销售金额的比例系数；vi：第i个月月初的库存量；si：第i个月的销售量；5、模型的建立与求解5.1问题一的模型建立与求解（1）模型一的建立1）建立目标函数通过对问题一的深入分析，我们总结出问题一是一个最优化问题，即在满足该厂生产能力不受限制、该厂允许期货销售和该厂能够在上半年把生产的制砖机全部销售完三个条件下，最终使得该厂的利润最大。根据“利润=销售额-成本”的公式，结合题目所给的信息可知，该厂上半年的销售额等于制砖机的销售价格乘以上半年总的销售量，即p0xi6i=1(5-1)该厂上半年的成本包括生产成本（生产成本又包括固定成本和可变成本）、销售费用和贮存费，即a0+0.5*xi2*(1+ai)+k*p0*xi+0.1*vi666i=1i=1i=1因此，目标函数可以表示为：(5-2)p0xi-(a0+0.5*xi2*(1+ai)+k*p0*xi+0.1*vi)Max6666i=1i=1i=1i=1(5-3)2）建立约束条件根据题目所给的信息，结合实际情况可知，第i+1个月月初的库存量等于第i个月月初的库存量加上第i个月生产的数量减去第i个月的销售数量，即vi+1=xi+vi-si(5-4)此外，1月份月初的库存量v1=49，并且每个月的生产计划xi0。3）建立非线性优化模型根据题目可知，该厂上半年的销售额、固定成本和销售成本都是常数，因此要使得该厂的利润最大，也就是使得该厂的成本最低，所以目标函数可以简化为(0.5*xMin6i=12i*(1+ai)+0.1*vi)(5-5)综上所述，为使得该厂上半年的利润最大，即成本最低，制定该厂上半年的生产计划的非线性优化模型如下：(0.5*x目标函数：Min6i=12i*(1+ai)+0.1*vi)约束条件：s.tvi+1=xi+vi-si；v1=49；xi0；（2）模型一的转化通过查阅相关的文献资料，结合题目的背景和模型一的特点，我们认为该厂的生产计划问题属于多阶段决策问题，因此该非线性优化模型可以通过转化为动态规划模型，然后通过MATLAB软件编程来进行求解。1）动态规划的基本模型在实际应用中，要构造一个标准的动态规划模型，通常需要采用以下几个步骤：l划分阶段：按照问题的时间或空间特征，把问题分为若干个阶段，这些阶段必须是有序的或者是可排序的(即无后向性)，否则，应用无效。l选择状态：将问题发展到各个阶段时所处的各种客观情况用不同的状态表示，称为状态。状态的选择要满足无后效性和可知性，即状态不仅依赖于状态的转移规律，还依赖于允许决策集合和指标函数结构。l确定决策变量与状态转移方程：当过程处于某一阶段的某个状态时，可以做出不同的决策，描述决策的变量称为决策变量。在决策过程中，由一个状态到另一个状态的演变过程称为状态转移。l写出动态规划的基本方程：动态规划的基本方程一般根据实际问题可分为两种形式，逆序形式和顺序形式。动态规划基本方程的逆序形式为：Dk(v,fk(xk)=ukoptxk)k(xk,uk)+fk+1(xk+1)k=n,n-1,L,2,1件件件件件fn+1(xn+1)=0件fn(xn)=vn(xn,un)其中第k阶段的状态为xk，其决策变量uk表示状态处于xk+1的决策,状态转移方程为xk+1=Tk(xk,uk)，k阶段的允许决策集合记为Dk(xk)，vk(xk,uk)为指标函数。当求解时，由边界条件从k=n开始，由后向前逆推，逐阶段求出最优决策和过程的最优值，直到最后求出f1(x1)，即得到问题的最优解。类似地，动态规划基本方程的顺序形式为：fk(xk+1)=件件件件件v,optk(xk+1,uk)+fk-1(xk)k=1,2,L,n-1,nukDk(xk+1)f0(x1)=02）建立动态规划模型=0,4971vv综上，建立该厂的生产计划的动态规划模型的步骤如下：l划分阶段：将该厂上半年的生产计划按照时间进行划分，一个月份为一个阶段，因此一共是6个阶段。l选择状态变量：将该厂每个月月初的库存量vi设为状态变量。l确定决策变量与状态转移方程：将该厂每个月生产制砖机的数量xi设为决策变量；同时根据题意可得状态转移方程为vi+1=xi+vi-si。l写出动态规划的基本方程：本文采用的是逆序法构建动态规划模型，i其中指标方程为0.1*vi+0.5*x2*(1+ai)，表示第i个月的贮存费和变动成本；因此，动态规划的基本方程如下：0.1*vi+0.5*xi2*(1+ai)+fi+1(xi+1),ifi(xi)=min=6,L,2,1件件件件件综上，最终的转化成的动态规划模型如下：0.1*vi+0.5*xi2*(1+ai)+fi+1(xi+1),ifi(xi)=min=6,L,2,1vi+1=xi+vi-siv1=49v7=0xi03）基于MATLAB的动态规划模型求解对于构建的动态规划模型，若用人工进行计算，则非常复杂，又很繁琐。如果基于MATLAB的动态规划逆序算法，同时通过MATLAB软件编程进行求解，就可以大大提高计算速度和准确率。其中，动态规划逆序算法的MATLAB程序详见附录8.1。根据上述建立的动态规划模型，包括阶段指标函数、状态转移方程和基本方程，本文写出下面3个M函数文件，以便在求解时可以进行调用。ldecisfun.m文件%decisfun.m阶段状态变量决策变量阶段指标14933603.8524034639.834231580.8表5.1问题1计算结果functionu=decisfun(k,s,u)%在阶段kk由状态变量x的值求出其相应的决策变量所有的取值b=42,32,41,67,25,29;ifk=6u=b(k)-s;%6月初的库存量和6月份生产的数量等于6月份的销售量elseu=0:236;%u表示某个月的生产数量，所以不能取负值；%每个月的生产数量不超过上半年的销售总数量减去1月初的库存量，即u不大于236endltransfun.m文件%transfun.mfunctions_next=transfun(k,s,u)%状态转移函数b=42,32,41,67,25,29;%销售数量矩阵s_next=s+u-b(k);%状态转移方程lsubobjfun.m文件%subobjfun.mfunctionv=subobjfun(k,s,u)%阶段目标函数a=0.1,0.1,0.2,0.2,0.3,0.3;%建材价格增长矩阵ifs0%由于该厂可以期货销售，因此库存量可能取负值，当库存量取负值时，该厂的贮存费为零v=0.5*u2*(1+a(k)+0.1*s;elsev=0.5*u2*(1+a(k);end在MATLAB主窗口输入每一个阶段状态变量的所有可能取值，然后直接调用dynprog.m（详见附录8.1）进行计算，计算结果如表5.1所示（详见附录8.2）。月份123456生产量（台）33343131292943231579.85-429546.656029546.65由表5.1可知，在生产能力没有限制的前提下，该厂生产制砖机的月进度表即生产计划如表5.2所示。表5.2生产计划表15.2问题二的模型建立与求解（1）模型二的建立在问题一中，我们构建的模型是在假设该厂的生产能力没有受到限制的前提下，但是结合实际情况，任何工厂的生产能力都是受到限制的，比如说原材料的短缺等等，因此问题二就是在问题一的基础上增加了一个约束条件，限制该厂每月的生产能力不大于33台，即：xi33(5-6)因此，在限制该厂生产能力的条件下，该厂生产计划的非线性优化模型如下：(0.5*x目标函数：Min6i=12i*(1+ai)+0.1*vi)约束条件：s.tvi+1=xi+vi-si；v1=49；0xi33；同样，可以将该非线性优化模型转化成动态规划模型，转化的步骤在前面的文本中已给出，因此这边只给出最终转化后的动态规划模型，如下：0.1*vi+0.5*xi2*(1+ai)+fi+1(xi+1),ifi(xi)=min=6,L,2,1阶段状态变量决策变量阶段指标14933603.8524033602.9534131580.743132617.55-429546.65表5.3问题2计算结果vi+1=xi+vi-siv1=49v7=00xi30（2）模型二的求解由于该动态规划模型只是增加了对决策变量的一个范围约束，因此阶段指标函数、状态转移方程和基本方程都没有发生改变，只需要改变decisfun.m文件中u的取值范围即可。改变后的decisfun.m（字体加粗且有下划线的地方即为改变的地方）为：%decisfun.mfunctionu=decisfun(k,s,u)%在阶段k由状态变量x的值求出其相应的决策变量所有的取值b=42,32,41,67,25,29;ifk=6u=b(k)-s;%6月初的库存量和6月份生产的数量等于6月份的销售量elseu=0:33;%u表示某个月的生产数量，所以不能取负值；该厂每个月的生产能力限于33台end同样，在MATLAB主窗口输入每一个阶段状态变量的所有可能取值，然后直接调用dynprog.m（详见附录8.1）进行计算，计算结果如表5.3所示（详见附录8.3）。月份123456销售量423241672529生产计划333431312929库存量49404232-40价格上涨幅度0.10.10.20.20.30.3月份123456生产量（台）3333313229296029546.65由表5.3可知，在生产能力限于33台的前提下，该厂生产制砖机的月进度表即生产计划如表5.4所示。表5.4生产计划表25.3结果分析（1）问题一的结果分析根据求得的生产计划表，如表5.1所示，结合销售量、库存量和价格上涨的具体信息，如表5.5所示，我们对求得的结果根据变量的变化趋势进行了分析，如图5.1所示。表5.5变量的具体数据1807060504030201000.350.30.250.20.150.10.05销售量生产计划1库存量价格上涨-101234560图5.1变量变化的趋势1（注：图中价格上涨的数据是以次坐标轴为单位，其余都是以主坐标轴为单位）月份123456销售量423241672529生产计划333331322929库存量49404131-40价格上涨幅度0.10.10.20.20.30.3根据图5.1中变量变化的趋势可知，该厂上半年的销售数量的波动比较大，而生产计划比较稳定，但是总体呈现下降的趋势。这是由于建材成本的价格上涨幅度一直处于增长的状态，因此为了节约生产成本，该厂应该尽可能多的在建材成本较低的时候进行生产，这符合实际情况。由于库存量受到销售量和生产计划的共同作用，因此库存量的波动也比较大，特别是在5月份，出现了库存量小于零的状况，这说此时该厂进行了期货销售。（2）问题二的结果分析在该厂的生产能力受到限制后，同样根据求得的生产计划表，如表5.3所示，结合销售量、库存量和价格上涨的具体信息，如表5.6所示，我们对求得的结果根据变量的变化趋势进行了分析，如图5.2所示。表5.6变量的具体数据2807060504030201000.350.30.250.20.150.10.05销售量生产计划2库存量价格上涨-101234560图5.2变量变化的趋势2（注：图中价格上涨的数据是以次坐标轴为单位，其余都是以主坐标轴为单位）图5.2中变量变化的趋势和图5.1中变量变化的趋势基本一致，有所不同的是图5.2求出来的结果是在该厂的生产能力受到限制后求得的，因此该厂的生产计划不能超过每月33台。从图5.2中的可以看出，该厂在第1、2月份生产的制砖机数量达到了最大生产能力33台，主要是因为1、2月份的建材成本最低。由于建材成本呈现上升的趋势，因此将更多的制砖机在建材成本较低的时候进行生产，就可以降低变动成本，最终提高该厂的利润，因此求得的结果是符合实际的。6、模型的评价与改进6.1模型的优点动态规划是解决多阶段决策过程最优化问题的一种方法，本文针对制砖机生产计划多阶段决策问题的特点，把多阶段过程转化为一系列单阶段问题，利用各阶段之间的关系，逐个求解，提出了解决问题的最优化原理。动态规划模型可以通过MATLAB编程得到最终最优结果，避免了人工求解时的大量计算，节省时间且提高计算准确率。6.2模型的缺点在模型中，没有考虑工厂的实际库存，但是在实际生产中，每个工厂在每阶段的生产中都存在最大库存量；模型是严格基于2011年上半年的销售预测的基础上建立得到的，没有一定的弹性，如果预测略有偏差或者每阶段的预测有波动，此模型在一定程度上就不适用了；在模型中，假设只在月末进行销售，但在实际生产中，每天都存在商品交易6.3改进的方向针对上述模型所存在的缺点，后期的改进方向可定为三个方面：（1）在模型中增加一个变量每阶段该厂分配给制砖机的库存量；（2）将销售预测改为变动值，用预测模型代替销售预测值；（3）将库存费用以天为单位进行计算，使构建的模型更合理。7、参考文献1褚洪生，杜增吉等.MATLAB7.2优化设计实例指导教程.机械工业出版社，2007.2孙晓君.基于MATLAB的动态规划逆序算法的实现J.纺织高校基础科学学报，2002,15（1），38-41.3史历，夏先进.基于MATLAB解决生产库存问题J.制造业自动化，2010,32（10），51-52.4于斌,刘姝丽,韩中庚.动态规划求解方法的MATLAB实现及应用J.信息工程大学学报，2005,6（3），95-98.5丰玮，吴凤平.对生产计划问题DP模型的进一步探讨J.河海大学学报，1998,26（3），46-49.6孙宝，王希云.基于MATLAB的动态规划常用算法的实现J.太原师范大学学报，2008,7（4），26-30.8、附录8.1动态规划逆序算法的MATLAB程序functionp_opt,fval=dynprog(x,DecisFun,SubObjFun,TransFun,ObjFun)k=length(x(1,:);x_isnan=isnan(x);t_vubm=inf*ones(size(x);f_opt=nan*ones(size(x);d_opt=f_opt;tmp1=find(x_isnan(:,k);tmp2=length(tmp1);fori=1:tmp2u=feval(DecisFun,k,x(i,k);tmp3=length(u);forj=1:tmp3tmp=feval(SubObjFun,k,x(tmp1(i),k),u(j);iftmp=t_vubm(i,k)f_opt(tmp1(i),k)=tmp;d_opt(tmp1(i),k)=u(j);t_vubm(i,k)=tmp;endendendforii=k-1:-1:1tmp10=find(x_isnan(:,ii);tmp20=length(tmp10);fori=1: