公司收益问题

1、公司收益问题摘要本文运用了规划模型解决托运公司对客户申请量的的批复计划，使得托运公司能最大限度的获得高收益。在问题一中建立规划模型对托运公司在一天中的托运计划作出优化安排，在A类6500kg，B类5000kg，C类4000kg，D类3000 kg的客户申请量的前提下，运用lingo软件求解该规划问题，得出最大获利为：40232元。具体的决策如下：(a) : 3辆卡车中，用一辆卡车只托运 C：服装类 kg，而A：鲜活类、B：禽苗类、D：其他类货物的托运量为都为0 kg；(b) : 用一辆卡车托运 A：鲜活类 kg，托运 C：服装类 kg，其余两类货物的托运量为0 kg；(c) : 用一辆卡车托运

2、 B：禽苗类 kg，托运 C：服装类 kg，其余两类货物托运量为0 kg。问题二中，因公司以往的客户托运申请量是随机的，数据的波动性较大，在预测分析时, 经常要碰到趋势随机型数据序列(既具有一定的趋势性又具有较大随机性的数据序列) 的预测问题。考虑到用一般的时间序列方法对公司未来7天的申请量进行时，由于数据的不平稳性及随机性，预测结果会出现较大的误差，在本文中采用了，随机型时间序列模型进行预测。随机型时间序列模型的思想是：将随机型时间序列看成由发展趋势曲线和围绕趋势曲线的波动二部分迭加而成, 用线性回归等方法确定数据的发展趋势曲线，围绕趋势曲线划分成若干个状态。用马尔可夫概率矩阵预测分析状态的

3、转移规律，确定数据的可能转移状态，得到预测值。最后预测得出未来7天客户的托运申请量如下表所示：天数A：鲜活类B：禽苗类C：服装类D：其他类12874336878003464218653531677435223223235346364347742074354860953468521343548596134616211235495874346072120354958283549问题三中，是与问题一类似最优规划问题，最后求解得未来7天，各天的最大收益为：天数1234567各天的最大收益（元）40878407414063640554405154048940475【关键词】规划模型 lingo软件 随机

4、型时间序列预测模型 马尔可夫概率矩阵一 问题重述某货运公司拥有3辆卡车，每辆卡车载重量均为8000kg可载体积为9.084m3，该公司为客户从甲地托运货物到乙地，收取一定费用。托运货物可分为四类，A：鲜活类，B：禽苗类，C：服装类，D：其他类，公司有技术实现四类货物任意混装。平均每类每kg所占体积和相应托运单价如下表：类别A：鲜活类B：禽苗类C：服装类D：其他类体积(m3/kg)0.00120.00150.0030.0008托运单价(元/kg)1.72.254.51.12托运手续是客户首先向公司提出托运申请，公司给予批复，客户根据批复量交货给公司托运。申请量与批复量均以公斤为单位，例如客户申请

5、量为1000kg，批复量可以为01000kg内的任意整数，若取0则表示拒绝客户的申请。问题一、如果某天客户申请量为：A类6500kg，B类5000kg，C类4000kg，D类3000kg，如果要求C类货物占用的体积不能超过B、D两类体积之和的三倍（仅在问题一中作此要求）。公司应如何批复，才能使得公司获利最大？问题二、每天各类货物的申请总量是随机量，为了获取更大收益，需要对将来的申请总量进行预测。现有一个月的数据（见附件一），请预测其后7天内，每天各类货物的申请量大约是多少？问题三、一般，客户的申请是在一周前随机出现的，各类申请单位立即批复，批复后即不能更改，并且不能将拒绝量（即申请量减批复量）

6、累计到以后的申请量。根据对下周7天中各类货物申请量的预测，估计这7天的收益各为多少？附件一 某月申请量数据表（单位：kg）日期A类B类C类D类总计11601284549262239116112542128332871243113683189044884447275013575444394554299614841347351703292850884378140976323234972829359313151737622613893211786478116769216706183716667918971391806417501310210373735083386593816641111807445

7、153171459130341216282636311277571513313172334714226244111861142584385445201373123311515513556349423651096616247926592918266010716171199433528603078114721841482882551436361618019244940842008308111622202026199958223204130512116902889284013188737223374217528934083125252320152510112138339479242480340916

8、631773932525850372927362519983426224934894552605016340271674317287944710183502836664568555211791496529202940151195323932039030123836669552257917035二 模型假设1货物在装车后由于挤压而产生的体积变化不大。2该公司的3辆卡车除载重量（8000kg）和容积（9.084m3）都相同外，卡车的车况也相同，可认为是3辆同样的车。3对3辆卡车分别标上记号，称为卡车1、卡车2、卡车3。4公司进行货物装车，在运用混装技术时对卡车的装载能力没有影响。5每天各类货物的申

9、请总量是随机量。6假设A类、B类、C类、D类货物分别称为第1类、第2类、第3类、第4类货物。7. 公司每次的批复量均为整数。8. 问题三中客户的申请在一周前随机出现，客户的申请立即得到批复且批复后不能再更改，不能将拒绝量累计到以后的申请量。三 全局符号说明以下符号定义中： ：公司拥有的卡车数量， ：每辆卡车的载重量， ：每辆卡车的容量， ：第类货物单位重量kg所占的体积。 ：第类货物从甲地托运到乙地单位重量kg的运价。 ：卡车装运第类货物的重量。 ： A类、B类、C类、D类货物的托运申请量。 ：公司对第类货物的拒绝量。四 问题分析1对问题一 ，要求的是是公司在客户申请量已知的情况下，应该如何批

10、复申请量才能使收益最大，我们通过建立优化模型，利用lingo软件进行求解从而给出最优方案。2对问题二，考虑到客户托运申请量的随机性，我们运用随机型时间序列预测方法，通过matlab编程可以得出未来一周的预测量。通过处理和分析问题中所提供的原始数据，发现在过去的某一30天内各天中，客户的各类货物的托运申请量是随机的，而且数据的波动性很大，由此，若想单靠运用一般线性回归、曲线的拟合及一般的时间序列分析法去进行预测，则最后得出的预测结果的误差将会变得相当大，为此，在以下的建模预测过程中，引进了随机型时间序列预测法，该方法在预测既具有一定的趋势性又具有较大随机性的数据时，精确度较高。该随机型时间序列预

11、测模型的主要思路及步骤是： 预测趋势曲线的确定。 状态区间的划分。 状态转移概率的确定。预测值计算。3对问题三，要求未来一周各天的收益，实际是结合问题二的数据和题目的要求，对公司在未来7天内的收益进行规划和预测，运用与问题一类似的方法，再一次建立优化模型，通过lingo运行得到结果。五 模型建立与求解1.问题一1.1背景：已知客户在一天中提出的托运申请量为：A类6500 kg，B类5000kg，C类4000kg，D类3000kg。公司有技术实现四类货物任意混装由假设1，货物在装车后体积的变化不大，即可认为装车后货物的体积不变。公司应如何对客户提交的货物申请量进行批复才能使得公司的获利最大。现将

12、A、B、C、D四类货物的单位重量kg所占的体积、相应托运单价、托运申请量列表：类别 A：鲜活类B：禽苗类C：服装类D：其他类体积(m3/kg)：托运单价元/kg元/kg元/kg元/kg托运申请量kgkgkgkg1.2 问题一模型建立与求解为使问题简化，对3辆卡车分别标上记号，称为卡车1、卡车2、卡车3，三辆卡车可认为是相同的车。公司在决定批复量的多少时，首先须要对各辆卡车分别装运多少重量的哪类货物作出决策，使得获利最大。1.21 ：公司的获利目标函数为：在决策时必须满足各种条件限制：1.22：首先是卡车载重量，卡车的容量的限制，则：公司对各类货物的批复量不能超过客户的托运申请量：A类kg，B类

13、kg，C类kg，D类kg，则： C类货物占用的体积不能超过B、D两类体积之和的三倍，则：由可得问题一模型：运用lingo8.0进行最优化求解，（具体程序及运行结果详见附件1）：求解结果：公司的最大获利为：40232元具体的决策如下表一：A：鲜活类B：禽苗类C：服装类D：其他类卡车1装运情况（单位：千克） kg kg kg kg卡车2装运情况（单位：千克） kg kg kg kg卡车3装运情况（单位：千克 kg kg kg kg由表一可得公司获利最大时的决策方案为：(a) : 3辆卡车中，用一辆卡车只托运 C：服装类 kg，而A：鲜活类、B：禽苗类、D：其他类货物的托运量为都为0 kg；(b)

14、: 用一辆卡车托运 A：鲜活类 kg，托运 C：服装类 kg，其余两类货物的托运量为0 kg；(c) : 用一辆卡车托运 B：禽苗类 kg，托运 C：服装类 kg，其余两类货物托运量为0 kg。由决策方案（表一），可得公司对客户的拒绝量如下表二：表二 单位：kg拒绝量：A：鲜活类B：禽苗类C：服装类D：其他类 kg kg kg kg表二说明：公司对 A：鲜活类货物的拒绝量为40 kg；对 B：禽苗类货物和 C：服装类货物的拒绝量为0 kg，即对这两类货物托运量（批复量）等于客户申请量；对 D：其他类货物的拒绝量为3000 kg，即批复量为0 kg，公司不托运该类货物。2. 问题二模型建立与求解

15、建立随机型时间序列模型对未来7天内各天中4类货物的申请量进行预测。建模步骤为：预测趋势曲线的确定。状态区间的划分。状态转移概率的确定。 预测值计算。在确定预测趋势曲线时，先大概的确定能反映总体数据走势的曲线，但此时的曲线并不是作为最终的预测用。在趋势曲线确定的前提下，主要采用最小二乘法确定各边界线，划分状态区间。在状态区间划分好后，建立马尔可夫概率矩阵，以确定状态转移概率矩阵。最后，通过考察转移概率矩阵, 明确了系统未来各状态的转移概率后, 就可确定预测值, 即由各状态的中点值与相应各状态的概率乘积之和求得。2.1 预测趋势曲线的确定 设趋势曲线的一般化方程表示为：2.2 状态的划分以趋势曲线

16、为基准，划分为与趋势曲线相邻的若干个条形区域，每一个条形区域构成一个状态。以下建模过程中划分为个状态，每一个状态用表示：其中为时间的函数，可能是时间的函数，随时序变化，即状态具有动态性。状态数目和、的确定，可根据研究对象和观测数据数目来确定。对A类货物进行预测时，选择，即在以下分析A类货物的申请量数据时划分为4个状态。(a)上下极限边界线的确定上下极限边界线的确定原则是： 所有观测点（数据点）都必须落在上极限边界线以下和下极限边界线以上。 使上下极限边界线与趋势曲线对应各点之差的为最小。上下极限边界线的形状有两种情况，一是上下极限边界线与趋势曲线形状相同，则上下极限边界线可由趋势曲线平移而得。

17、一是上下极限边界线与趋势曲线的形状不同，则根据散点图的形式选择适应的线性或非线形模型。根据A类货物的申请量数据情况，划分为四个状态区间，设A类货物的上下极限边界线及划分状态区间时的上下边界线形状与趋势曲线的形状相同，则A类货物的各边界线可由其趋势曲线通过平移得到。 设上极限边界线为：由边界线曲线确定原则和可得、：其中为对应于时的原始数据值或观测点。由式可知，只要求得最小的值，便得最小的值，值同时必须满足式，即对所有的数据（观测数据或原始数据）都能使式成立。可得： 设下边界线为：同理，由边界线曲线确定原则和可得、式：同理，由、可得：(b)状态区间边界线的确定将趋势曲线以上各区间的边界线定义为上边

18、界线；将趋势曲线以下各区间的边界线定义为下边界线。上下边界线采用最小二乘法来确定，即用落入趋势曲线以上的观测数据对上边界线的残差平方和为最小来确定上边界线；用落入趋势曲线以下的观测数据对下边界线的残差平方和为最小来确定下边界线。若上下边界线形状与趋势曲线形状相同则可通过由趋势曲线平移而得到。在前面已设定A类货物的上下极限边界线及划分状态区间时的上下边界线形状与趋势曲线的形状相同，则A类货物的各边界线可由其趋势曲线通过平移得到。设上边界线为：则由最小二乘法原理，得：式中为落入趋势曲线以上的各观测值。根据微分方程求极值的方法，应满足：则可得：其中P + 1 为落入趋势曲线以上的观测数据个数。设下边

19、界线为：则由最小二乘法原理，得：根据微分方程求极值的方法，应满足：则可得：式中q+ 1 为落入趋势曲线以下的观测数据个数。(c)区间状态的确定确定了上下极限边界线和上下边界线后, 则各区间的上、下边界也就确定了。若划分四个区间, 则分别为：因为上下边界线和上下极限边界线曲线的形状都与趋势曲线形状相同, 则所划分的四个区间分别为:(d)状态转移概率的确定建立马尔可夫概率矩阵。对阶马尔可夫序列, 从状态经过步转移到状态的概率由下式表示：式中 ：表示由状态经过m 步转移到状态的概率(本文中一步代表天年) ; ：为由状态经过步转移到状态的观测数据样本数; ：为处于状态的样本数。则状态转移概率矩阵为:状

20、态转移概率反映了系统内部各状态之间转移的规律。考察状态转移概率矩阵，可预测系统未来状态的可能转向。考察一步转移概率矩阵，设预测对象处于，则考察矩阵的第Q行，第Q行的各概率反映了转台转向其他各状态的概率。2.4预测值的计算通过考察转移概率矩阵, 明确了系统未来各状态的转移概率后, 就可确定预测值, 即由各状态的中点值与相应各状态的概率乘积之和求得:2.5 模型求解现将托运公司某月的申请量列表（具体数据见问题重述）。月申请量数据表（单位：kg）日期A类B类C类D类总计116012845492622391161125421283328712431136829202940151195323932039

21、030123836669552257917035预测趋势曲线主要反映观测数据总体的变化规律，可以用线性回归、非线性回归等方法来确定预测曲线。现用matlab软件分别作出上表4类货物的散点图，进行拟合后确定各类货物的预测趋势曲线。在以下建模中，经拟合后，采用一次线性曲线做为各类货物的趋势曲线：在matlab中分别作出A、B、C、D四类货物的申请量数据散点及折线图，在图形窗口的【tools】菜单中用linear选项进行拟合，可得四类货物在30天内一次线性趋势曲线（直线（linear）的图表1、2、3、4按顺序从上到下如下：图表 1、2、3、4图表 2、2、3、4日期：天数由图表1、2、3、4，可看

22、出A、B、C、D四类货物的粗糙的趋势曲线，即，四类货物的趋势曲线分别为：根据式，即：由上5个式子，用MATLAB编程可得：首先对A类货物进行分析：由：可得，A类货物的状态区间：可得A类货物的状态图如下所示：根据状态图可知：落入，、状态的观测数据样本点数分别为、,和，由状态一步转移到状态，、的观测数据样本数分别为,和同理, 可计算的值。 可得各状态间的一步转移概率矩阵如下由参考文献可知在此模型中转移概率矩阵即，在此，不转移概率矩阵是一步转移概率矩阵的次方，由此可得：到步转移概率矩阵，再此，步转移概率矩阵见附件。由式对A类货物在未来7天内的申请量进行预测可得，A类货物的预测值为：天数1234567

23、预测申请量（kg）2874186522322074213421122120注：程序及结果在附件。现作出B、C、D类货物的状态图形如下： 同理，可得B、C、D 类货物申请量数据的各个状态及概率各阶概率转移矩阵，注：所编程序及结果见附件。最终求解得预测值为： 天数A：鲜活类B：禽苗类C：服装类D：其他类128743368780034642186535316774352232232353463643477420743548609534685213435485961346162112354958743460721203549582835493问题三模型建立与求解3.1问题三符号说明：以下： ：第类货物

24、单位重量kg所占的体积。 ：第类货物从甲地托运到乙地单位重量kg的运价。 ：每辆卡车的载重量， ：每辆卡车的容量， ：卡车在一天中装运第类货物的重量。 : 由前面的随机型时间序列模型对未来7天中第类货物在第天的预测的申请量。3.2 问题的分析问题三中的要求是由问题二中预测出的7天申请量的基础上分别估计各天的收益为多少，在问题一中，其中的一个约束条件要求C类货物占用的体积不能超过B、D两类体积之和的三倍，在问题三中并未作要求，因此，可对问题一中建立的规划模型去掉次约束条件，可得问题三的模型。由前面的随机型时间序列模型对原有数据进行处理分析后，预测得到在未来7天内各天4类货物的申请量，如下表所示：

25、天数A：鲜活类（）B：禽苗类（）C：服装类（）D：其他类（）12874336878003464218653531677435223223235346364347742074354860953468521343548596134616211235495874346072120354958283549单位：kg在公司卡车数量为以及卡车的载重量不变的情况下，估算未来七天各天的收益。公司未来七天中每天的收益可以用如下的目标函数表示：利用LINGO软件对此模型分别对7天的预测数据进行求解可得各天的最大收益为：天数1234567各天的最大收益（元）40878407414063640554405154048

26、940475具体决策安排见附录件。六 模型评价与推广1 模型的优点将既有趋势性又有较大波动性的时间数据序列看成由发展趋势曲线和围绕趋势曲线的波动二部分迭加而成, 建立了较为成熟的数学模型，可信度高 。所建模型是通过对大量数据统计进行分析，可解决类似问题，有较强的推广性。用Matlab软件进行编程，从而求出马尔可夫概率矩阵，减少了人眼观察和手工循环整合的误差，使整个模型更加科学，评价的结果更加可信。实用性强，便于具体操作。2模型的缺点随机型时间序列的预测趋势曲线只取线性规划的一 次曲线，预测值能反映未来申请量的平均情况，但没有经过和其他预测方法进行比较分析。由于时间的关系，只考虑到上下极限边界线

27、及上下边界线平行和预测趋势曲线形状相同，没有分析形状不同的情形。3随机型时间序列模型的推广从问题的预测模型分析，所得到的预测数据比实际数据平稳程度要高出很多，用随机型时间序列模型进行预测，在该模型中，其更能反映的是未来7天之内各天的申请量的平均水平，而实际的数据会比预测值的波动性要大。但在本文中，用该模型对未来时期内的申请量进行预测，使公司能较好的对未来时间段的平均收益作出较为准确的预测，以使公司能对未来的营销作出规划。该模型不仅适用于公司收益问题，还可以推广到其他很多诸如此类的优化问题，如人力资源分配等问题。七 问题延伸本文的模型是在运输条件大为简化的情况下建立的，当问题中给出更多的条件限制

28、时，如考虑到运费和车辆定期维护所需的时间、车辆不仅3辆，而是多辆卡车、公司或客户对托运问题存在更多的考虑时，所建立的模型应变得复杂而更符合实际。此时模型要提出修改。八 考文献1 姜启源，谢金星，叶俊，数学模型（第三版），北京，高等教育出版社，82130页，333338页，2003.8.2 盛骤，谢式千，潘承毅，概率论与数理统计（第三版），北京，高等教育出版社，354367页，2001.12.3 王沫然，MATLAB与科学计算（第二版），北京，电子工业出版社，2003.9.4 刘卫国，陈昭平，张颖，MATLAB程序设计与应用，北京，高等教育出版社，2002.6.5 魏巍，应用数学工具箱技术手册，

29、北京，国防工业出版社，2004.1.九 附录注：以下是其中一部分主要程序，其他程序及运行结果在附件。1求解问题一的程序及运行结果程序modle:sets:class/1.4/:b,v;car/1.3/:;links(car,class):x,c;endsets data:v=0.0012 0.0015 0.003 0.0008;c=1.7 2.25 4.5 1.12 1.7 2.25 4.5 1.12 1.7 2.25 4.5 1.12;b=6500 5000 4000 3000;enddatamax=sum(links:c*x);for(car(i):sum(class(j):x(i,j)&

30、lt;=8000);for(car(i):sum(class(j):x(i,j)*v(j)<=9.084);sum(car(i):x(i,3)*v(3)-3*(x(i,2)*v(2)+x(i,4)*v(4)<=0;for(class(j):sum(car(i):x(i,j)<=b(j);for(links:gin(x);运行结果：Global optimal solution found at iteration: 7 Objective value: 40232.00 Variable Value Reduced Cost B( 1) 6500.000 0.000000 B

31、( 2) 5000.000 0.000000 B( 3) 4000.000 0.000000 B( 4) 3000.000 0.000000 V( 1) 0.1200000E-02 0.000000 V( 2) 0.1500000E-02 0.000000 V( 3) 0.3000000E-02 0.000000 V( 4) 0.8000000E-03 0.000000 X( 1, 1) 0.000000 -1.700000 X( 1, 2) 0.000000 -2.250000 X( 1, 3) 3028.000 -4.500000 X( 1, 4) 0.000000 -1.120000 X

32、( 2, 1) 6460.000 -1.700000 X( 2, 2) 0.000000 -2.250000 X( 2, 3) 444.0000 -4.500000 X( 2, 4) 0.000000 -1.120000 X( 3, 1) 0.000000 -1.700000 X( 3, 2) 5000.000 -2.250000 X( 3, 3) 528.0000 -4.500000 X( 3, 4) 0.000000 -1.120000 C( 1, 1) 1.700000 0.000000 C( 1, 2) 2.250000 0.000000 C( 1, 3) 4.500000 0.000

33、000 C( 1, 4) 1.120000 0.000000 C( 2, 1) 1.700000 0.000000 C( 2, 2) 2.250000 0.000000 C( 2, 3) 4.500000 0.000000 C( 2, 4) 1.120000 0.000000 C( 3, 1) 1.700000 0.000000 C( 3, 2) 2.250000 0.000000 C( 3, 3) 4.500000 0.000000 C( 3, 4) 1.120000 0.000000 Row Slack or Surplus Dual Price 1 40232.00 1.000000 2

34、 4972.000 0.000000 3 1096.000 0.000000 4 2472.000 0.000000 5 0.000000 0.000000 6 0.000000 0.000000 7 0.000000 0.000000 8 10.50000 0.000000 9 40.00000 0.000000 10 0.000000 0.000000 11 0.000000 0.000000 12 3000.000 0.000000问题二 2 程序及运行结果cleara=1601 2845 4926 22395421 2833 2871 2431890 4488 4447 2750443

35、9 4554 2996 14841703 2928 5088 43783232 3497 2829 3593376 2261 3893 21171167 6921 6706 18731897 1391 8064 17503737 3580 3386 59381807 4451 5317 14591628 2636 3112 77571723 3471 4226 24412584 3854 4520 13731551 3556 3494 23652479 2659 2918 26601199 4335 2860 30784148 2882 5514 36362449 4084 2008 3081

36、2026 1999 5822 32041690 2889 2840 13183374 2175 2893 40832015 2510 1121 38332480 3409 1663 1773850 3729 2736 25192249 3489 4552 60501674 3172 8794 47103666 4568 5552 11792029 4015 11953 23931238 3666 9552 2579;%a是数据矩阵xa=length(a);%a的列长度x=1:xa;m=1;%m是多项式拟合的次数b1=polyfit(x',a(:,1),m);b2=polyfit(x&#

37、39;,a(:,2),m);b3=polyfit(x',a(:,3),m);b4=polyfit(x',a(:,4),m);b(1,:)=b1;b(2,:)=b2;b(3,:)=b3;b(4,:)=b4;%统一为bclear b1 b2 b3 b4for i=1:4 y(i,:)=zeros(1,xa); for j=1:m+1 temp=b(i,j)*x.(m-j+1); y(i,:)=y(i,:)+temp; endendclear tempC=max(a-y');D=max(y'-a);p=sum(a-y')>=0);q=sum(y'

38、-a)>=0);A=sum(a-y').*(a-y')>=0)./p;B=sum(y'-a).*(y'-a)>=0)./q;for i=1:4 yhi(1,:)=y(i,:)-D(i); yhi(2,:)=y(i,:)-B(i); yhi(3,:)=y(i,:); yhi(4,:)=y(i,:)+A(i); yhi(5,:)=y(i,:)+C(i); subplot(2,2,i); plot(x,yhi,x,a(:,i); thf(1,i)=(yhi(1,xa)+yhi(2,xa)/2; thf(2,i)=(yhi(2,xa)+yhi(3,xa

39、)/2; thf(3,i)=(yhi(3,xa)+yhi(4,xa)/2; thf(4,i)=(yhi(4,xa)+yhi(5,xa)/2;endP1,1=1/6 1/6 2/6 2/6 0 3/10 3/10 4/10 3/7 4/7 0 0 1/6 3/6 2/6 0;P1,2=1/5 1/5 2/5 1/5 1/8 1/8 3/8 3/8 2/10 2/10 5/10 1/10 1/6 3/6 1/6 1/6;P1,3=5/8 1/8 2/8 0 1/9 2/9 4/9 2/9 2/8 5/8 0 1/8 0 1/4 1/4 2/4;P1,4=0 3/7 2/7 2/7 2/9 4/9

40、1/9 2/9 3/8 1/8 4/8 0 2/5 1/5 1/5 1/5;for i=1:4 for j=2:7 Pj,i=Pj-1,i*P1,i; endendbp=2,3,4,2;for i=1:4 for j=1:7 YC(j,i)=Pj,i(bp(i),:)*thf(:,i); endendclear i j bp xa程序及运行结果clear40153666;119539552;xa=length(a);x=1:xa;m=1;b1=polyfit(x',a(:,1),m);b2=polyfit(x',a(:,2),m);b3=polyfit(x',a(:,3

41、),m);b4=polyfit(x',a(:,4),m);b(1,:)=b1;b(2,:)=b2;b(3,:)=b3;b(4,:)=b4;for i=1:4 y(i,:)=zeros(1,xa); for j=1:m+1 temp=b(i,j)*x.(m-j+1); y(i,:)=y(i,:)+temp; endendclear tempC=max(a-y');D=max(y'-a);p=sum(a-y')>=0);q=sum(y'-a)>=0);A=sum(a-y').*(a-y')>0)./p;B=sum(y'

42、;-a).*(y'-a)>0)./q;for i=1:4 yhi(1,:)=y(i,:)-D(i); yhi(2,:)=y(i,:)-B(i); yhi(3,:)=y(i,:); yhi(4,:)=y(i,:)+A(i); yhi(5,:)=y(i,:)+C(i); subplot(2,2,i); plot(x,yhi,x,a(:,i); thf(1,i)=(yhi(1,xa)+yhi(2,xa)/2; thf(2,i)=(yhi(2,xa)+yhi(3,xa)/2; thf(3,i)=(yhi(3,xa)+yhi(4,xa)/2; thf(4,i)=(yhi(4,xa)+yhi

43、(5,xa)/2;endP1,1=1/6 1/6 2/6 2/6 0 3/10 3/10 4/10 3/7 4/7 0 0 1/6 3/6 2/6 0;P1,2=1/5 1/5 2/5 1/51/8 1/8 3/8 3/82/10 2/10 5/10 1/101/6 3/6 1/6 1/6;P1,3=5/8 1/8 2/8 01/9 2/9 4/9 2/92/8 5/8 0 1/80 1/4 1/4 2/4;P1,4=0 3/7 2/7 2/72/9 4/9 1/9 2/93/8 1/8 4/8 02/5 1/5 1/5 1/5;for i=1:4 for j=2:7 Pj,i=Pj-1,i*P1,i; endendbp=2,3,4