电力市场的输电阻塞模型 带matlab程序

1、电力市场的输电阻塞模型 摘要: 本文针对第一问，建立了多元线性均值回归模型，用于描述各线路上的有功潮流与各发电机组出力的关系。模型中利用均值定理，运用matlab软件求得解，此模型的误差较小，较符合实际情况。 对于第二问，我们引入基本补偿价概念，利用限上和限下的方法设计阻塞费用计算规则，在阻塞费用计算规则里面我们认为序外容量在低于对应报价的清算价上出力的时候，网方已经按清算价给予了发电商一部分的费用，这部分费用不作为经济补偿，即其是不包含在阻塞费用中，对于序外容量的经济补偿我们就用对应段容量的报价和清算价之差作为补偿价格来计算阻塞费用。 对于第三问，在考虑爬坡速率和负荷需求的前提下，根据电力市

2、场规则，考虑最小购电成本为目标函数，利用最优化的方法，我们得到了各机组的出力分配预案（如下表所示）。机组12345678出力1507918099.512514095113.9此时各机组所对应的潮流值（如下表所示）线路123456潮流值173.315796.8095-150.5051101.8151136.8380168.5319对于第四问，利用第一问求得的结果，检验第三问各机组的出力分配预案是否会引起输电阻塞，如果引起输电阻塞，根据输电阻塞管理原则并考虑经济原则，用计算机编程搜索调整预案，得出阻塞费用为405.8004元，新的方案见下表机组12345678出力1538822890.177515

3、29560.1116.1225此时各线路所对应的潮流值（如下表所示）线路123456潮流值165100.6408-154.8553105.5524131.5169159.5725 对于第五问，重复三，四问的工作，我们得出方案（如下表所示），阻塞费用为94.0377元。机组12345678出力1538822899.515215560.3117此时各线路对应的潮流值（如下表所示）线路123456潮流值173.461994.1143-154.9968104.9631135.3452160.4779问题的重述：电网公司在组织交易，调度和配送时,必须遵循电网“安全第一”的原则，同时制定一个电力市场交易规

4、则，按照购电费用最小的经济目标来运作。市场交易-调度中心根据负荷预报和交易规则制定满足电网安全运行的调度计划各发电机组的处理(发电功率)分配方案；在执行调度计划的过程中，还需实时调度承担AGC(自动发电控制)辅助服务的机组出力，以跟踪电网中实时变化的负荷。设某电网有若干台发电机组和若干条主要线路，每条线路上的有功潮流（输电功率和方向）取决于电网结构和各发电机组的出力。电网每条线路上的有功潮流的绝对值有一安全限值，限值还具有一定的相对安全裕度（即在应急情况下潮流绝对值可以超过限值的百分比的上限）。如果各机组出力分配方案使某条线路上的有功潮流的绝对值超出限值，称为输电阻塞。当发生输电阻塞时，需要研

5、究如何制订既安全又经济的调度计划。 电力市场交易规则：1. 以15分钟为一个时段组织交易，每台机组在当前时段开始时刻前给出下一个时段的报价。各机组将可用出力由低到高分成至多10段报价，每个段的长度称为段容量，每个段容量报一个价（称为段价），段价按段序数单调不减。在最低技术出力以下的报价一般为负值，表示愿意付费维持发电以避免停机带来更大的损失。2. 在当前时段内，市场交易-调度中心根据下一个时段的负荷预报,每台机组的报价、当前出力和出力改变速率，按段价从低到高选取各机组的段容量或其部分（见下面注释），直到它们之和等于预报的负荷，这时每个机组被选入的段容量或其部分之和形成该时段该机组的出力分配预案

6、（初始交易结果）。最后一个被选入的段价（最高段价）称为该时段的清算价，该时段全部机组的所有出力均按清算价结算。注释：（a） 每个时段的负荷预报和机组出力分配计划的参照时刻均为该时段结束时刻。（b） 机组当前出力是对机组在当前时段结束时刻实际出力的预测值。（c） 假设每台机组单位时间内能增加或减少的出力相同，该出力值称为该机组的爬坡速率。由于机组爬坡速率的约束，可能导致选取它的某个段容量的部分。（d） 为了使得各机组计划出力之和等于预报的负荷需求，清算价对应的段容量可能只选取部分。市场交易-调度中心在当前时段内要完成的具体操作过程如下：1、 监控当前时段各机组出力分配方案的执行，调度AGC辅助服

7、务，在此基础上给出各机组的当前出力值。2、 作出下一个时段的负荷需求预报。3、 根据电力市场交易规则得到下一个时段各机组出力分配预案。4、 计算当执行各机组出力分配预案时电网各主要线路上的有功潮流，判断是否会出现输电阻塞。如果不出现，接受各机组出力分配预案；否则，按照如下原则实施阻塞管理：输电阻塞管理原则：（1） 调整各机组出力分配方案使得输电阻塞消除。（2） 如果（1）做不到，还可以使用线路的安全裕度输电，以避免拉闸限电（强制减少负荷需求），但要使每条线路上潮流的绝对值超过限值的百分比尽量小。（3） 如果无论怎样分配机组出力都无法使每条线路上的潮流绝对值超过限值的百分比小于相对安全裕度，则必

8、须在用电侧拉闸限电。（4） 当改变根据电力市场交易规则得到的各机组出力分配预案时，一些通过竞价取得发电权的发电容量（称序内容量）不能出力；而一些在竞价中未取得发电权的发电容量（称序外容量）要在低于对应报价的清算价上出力。因此，发电商和网方将产生经济利益冲突。网方应该为因输电阻塞而不能执行初始交易结果付出代价，网方在结算时应该适当地给发电商以经济补偿，由此引起的费用称之为阻塞费用。网方在电网安全运行的保证下应当同时考虑尽量减少阻塞费用。现要解决的问题如下：（表格见附录1）1. 某电网有8台发电机组，6条主要线路，表1和表2中的方案0给出了各机组的当前出力和各线路上对应的有功潮流值，方案132给出

9、了围绕方案0的一些实验数据，试用这些数据确定各线路上有功潮流关于各发电机组出力的近似表达式。2. 设计一种简明、合理的阻塞费用计算规则，除考虑上述电力市场规则外，还需注意：在输电阻塞发生时公平地对待序内容量不能出力的部分和报价高于清算价的序外容量出力的部分。3. 假设下一个时段预报的负荷需求是982.4MW，表3、表4和表5分别给出了各机组的段容量、段价和爬坡速率的数据，试按照电力市场规则给出下一个时段各机组的出力分配预案。4. 按照表6给出的潮流限值，检查得到的出力分配预案是否会引起输电阻塞，并在发生输电阻塞时，根据安全且经济的原则，调整各机组出力分配方案，并给出与该方案相应的阻塞费用。5.

10、 假设下一个时段预报的负荷需求是1052.8MW，重复34的工作。问题的分析：1 初步分析题意，根据表1和表2所给的数据，我们利用数理统计知识和线性回归分析方法，拟建立多元线性均值回归模型，推算出各线路上有功潮流关于各发电机组出力的近似表达式；2 设计简明，合理的阻塞费用的计算规则，关键在于体现网方是怎样充分补偿发电厂的损失的，在这里，我们引入基本补偿价的概念，利用限上和限下的方法得出阻塞费用的计算规则；3 仅仅考虑电力市场交易规则和各机组的爬坡速率的影响，要得出下一时段各机组的出力分配预案，就是在各机组的总段容量不超出爬坡速率所限制的极限值的前提条件下，选取价格最低的全部机组的段容量总和，使

11、得这个总和等于负荷需求； 4 用第一问所得出来的有功潮流的表达式检查预案是否会引起输电阻塞，要是会发生输电阻塞，则根据输电阻塞管理原则并考虑电力市场规则的情况下调整预案得出一个新的方案，并根据设计的阻塞费用计算规则求出新方案对于预案的阻塞费用。模型的假设：1 各线路上的有功潮流只受到这个电网里面八个发电机组所发的电量的影响，与外界的电网无关，即电网结构不会改变；2 发电机组出力是按照爬坡速率一直增加或减少到符合出力分配方案的要求后，出力就是一直恒定不变的；3 发电机出力的增量是由低价格的容量增加到高价格的容量的，而出力的减量就是由高价格的容量减少到低价格的容量；符号的说明： ： 第i条线路的潮

12、流值（i1，2，6）； ： 序外容量要出力的段容量或其部分； ： 序内容量不能出力的段容量或其部分； ： 第i个机组的出力（i1，2，8）； ： 当前方案0中第i个机组的出力（i1，2，8）； ： 第i个机组在单位时间内的功率增量（i1，2，8）； ： 第i个机组的爬坡速率（i1，2，8）； ： 一个时间段（15分钟）； ： 对应段容量/对应段容量所在机组的爬坡速率； ： 对应段容量的前面容量总和变化到方案0的出力所需的时间；模型的建立：问题1 多元线性均值回归模型 根据题中所给统计数据，我们用回归分析方法，考虑线性分析模型，讨论得到：运用全回归分析虽然计算简单，但是最终的计算结果误差比较大。

13、所以我们建立多元线性均值回归模型。则有： （1）对（1）式的算法如下： （2） （3） 由表1和表2的132方案中知，为稀疏矩阵，元素只有一项i不为零（此项机组表现为相对0方案的变动），且同类的有4组，刚好8个机组的单独变动都有4组数，很明显我们可以均值就可求出i组的k的列向量，运用循环即可全部求得8组k的列向量，合并得k值。（程序附录2和3）系数矩阵k观察以上求出的系数k可以发现：k值比较小，而且相邻k值的变化也比较小。 求c时（程序见附录4），利用以上求得k值，有： （4） 与对应于132方案中的某一方案，由此运用循环当求36组，求均值可得。常数c所以写出y和p的近似表达式为：对上述方法进

14、行误差分析，由求出的表达式返回32组方案与原方案对比，得出y的相对误差图（如图一所示）： 如图一可知，此模型的误差较小。问题2对于第二问，考虑的是当阻塞发生的时候网方所要补偿给发电商的费用，这个费用称为阻塞费用S，一共包括两个部分：一部分是序内容量不能出力的部分，决定这部分阻塞费用的主要因素是不能出力的机组停发或少发电的度数以及少发一度电所需的补偿价格；另一部分是报价高于清算价的序外容量出力的部分，决定这部分阻塞费用的主要因素是序外容量机组多发电的度数以及多发一度电所需的补偿价格。模型一：限上和限下的方法根据上面的思想，我们应用限上和限下的方法1来计算阻塞费用，规则如下：1 当出现输电阻塞时，

15、在竞价中未取得发电权的发电容量（序外容量）按报价由低到高的顺序安排机组发电(限上),在竞价中本已经取得发电权的发电容量（序内容量）按报价由高到低的顺序安排机组停发或少发电(限下)。2 在竞价中未取得发电权的发电容量（序外容量）是在低于其对应报价的清算 价上出力的，这样就会使得发电商利益亏损而网方从中盈利，网方应该对这 些发电商给与经济补偿，这里补偿我们就用对应段容量的报价和清算价之差作为补偿价格；而那些在竞价中本已经取得发电权的发电容量（序内容量）却不能出力的发电商，由于不能出力或部分不能出力，使得这些发电商由于网方线路原因造成经济利益的损失，这时网方也应该给与这些发电商经济补偿，这里我们就以

16、清算价和其不能出力的对应的段容量的报价之差作为补偿价格。根据计算规则，我们可以列出阻塞费用S的式子： （5）其中 清算价对应段容量的报价对应段容量的报价清算价模型二：模型二是对于模型一所进行的改进，在模型一中可以发现：对于序内容量，当某个的段容量的报价等于它的清算价的时候，此时0，这表明这个段容量是得不到经济补偿的。由于每一个发电厂发的每度电都存在一个利润空间，这个利润是包含在报价里面的，按照模型一的结论表明发电厂按清算价作为报价的段容量的利润得不到经济补偿，显然这是不合理的。为了解决这个问题，我们在这里引入一个概念：基本补偿价B，根据文献资料 2 ，我们定 （6），于是我们就以清算价和其不能

17、出力的对应的段容量的报价之差再加上基本补偿价作为序内容量不能出力的发电商的经济补偿。于是式子（5）变为： （7）问题3由题意可知: （8）在下一个时间段的结束时刻各机组的出力大小为 （9）根据式子（8）（9），结合表1中方案0各机组的数据和表5中各机组的爬坡速率，相对于方案0，我们可以算出在下一个时间段的结束时刻各机组所允许能选取的出力范围（如下表31所示）表31 各机组所允许能选取的出力范围机组12345678最大出力1538822899.5152155102.1117最少出力875813260.5989560.163根据电力市场规则，我们对附录1中的表3，表4进行处理。（按照段价由低到高排

18、列如下表格群所示）表格群32 价格-800-610-607-590-560-505-500-500机组83651174段容量70110957570705055价格00000116120124机组12356571段容量00050151550价格146150152153159168170173机组54386146段容量054001001020价格180182183188189200203205机组72853426段容量520201501080价格210215233233245250251252机组15382571段容量0153001501530价格252253255258260283300302机

19、组68437824段容量1520100100610价格303305306308310312315318机组86735178段容量20101020100510价格320325330335348356360363机组24177321段容量21501034000价格380380396400405410415435机组64586234段容量20010150000价格489495500510520548800800机组12356784段容量408401010251按照以上表格，仅仅考虑电力市场规则，而不考虑各机组所允许的出力范围的时候，选取到第7个表格中的303价格位的段容量就能够满足负荷需求（这个价格

20、位之前的段容量全部选取），这时候我们得到各个机组的出力分别为：表33机组12345678出力1507918010012514095113.4由于受到爬坡速率的限制，第四机组的出力最大值为99.5，不能达到100，将第四组不能达到的0.5段容量按照电力市场规则选取还未选取的且对应段容量价格最低的，由此，我们得出下一个时段各机组的出力分配预案为（如下表34所示）：表34 下一时段各机组的出力分配预案（负荷需求为982.4MW时）机组12345678出力1507918099.512514095113.9问题4 将上题所得的各机组的出力分配预案代入第一问的模型二中求出的y和p的关系式中，得到各线路的潮

21、流值如下表41所示：表41 各线路的潮流值线路123456潮流值173.315796.8095-150.5051101.8151136.8380168.5319上表41与附录1中的表6相比较可以发现：线路1，5，6会引起输电阻塞。 图 二根据算法（如图二所示）编写程序（见附录6），由计算机搜索得到：改变后的出力分配方案如下表42所示：表42 各机组的出力分配方案机组12345678出力15388228901116.1225此时各线路的潮流值如下表43所示：表43 改变方案后各线路的潮流值线路123456潮流值165100.6408-154.8553105.552413

22、1.5169159.5725并得到对应的阻塞费用为405.8004元。问题5 利用第三问和第四问同样的方法，得出在下一个时段预报的负荷需求是1052.8MW时，各机组的出力分配预案如下（如下表51所示）表51 下一时段各机组的出力分配预案（负荷需求为1052.8MW时）机组12345678出力15081218.299.5135150102.1117此时各线路的潮流值如下表52所示表52 各线路的潮流值线路123456潮流值177.260894.6307-155.7527110.0612134.8521167.0870调整之后各机组的出力分配方案为：（如表53所示） 表53 调整后的各机组出力分

23、配方案机组12345678出力1538822899.515215560.3117各线路的潮流值为：（如表54所示） 表54 调整后的各线路的潮流值线路123456潮流值173.461994.1143-154.9968104.9631135.3452160.4779从表54可以知道，第1，5条线路都是超过限值的，但是都在安全裕度内，第一条线路超出限值的百分比为5.128，第五条线路超出限值的百分比为2.5342。并得出相应的阻塞费用为：94.0377元。模型的优缺点和改进方向：优缺点：本论文中的平均值分析法能在一定范围内很好地吻合数据，经误差分析亦可得相对误差0.0001左右，而且机组系数比较小

24、，能很好地说明此电网有比较好的稳定性，不会因一个机组的小变化而有大的改动，但是，此模型只适应机组出力都较大的情况下，在所有机组都比较小时，此模型有较大误差，但此范围并不在此题的计算范围内，此模型已经能非常很好地解决一般的电网调动。另外，本文中的计算较多应用matlab编程，对大量的数据计算提供很好的桥梁，针对每一个问题都有一个程序来解决，为以后的利用提供方便。桥梁改进方向：由于给出数据中有突变数据，所以在计算得出的相对误差中有6个是大于0.05的，当踢除这些值代以估算的均值，能更完美地完成数据的拟合。对于赔偿算法中有发电商的隐藏利润，此利润应该赔偿过加入更加合理准确，对此中的专业估计，是改进的

25、重点。参考文献：1 柯进,管霖 电力市场下的输电阻塞管理技术 电力系统自动化 第26卷 第14期 第四页 2002年2 西部水利水电投资项目介绍29/xbshj/xbbshd/xbbshd_indx.htm 访问时间：2004-9-19 韩忠诚 吕伯雄 金昌市招标项目 访问时间：2004-9-19附录1：表1 各机组出力方案 （单位：兆瓦，记作MW）方案机组123 4 5 6780120731808012512581.1901133.02731808012512581.1902129.63731808012512581.1903158.7773180801

26、2512581.1904145.32731808012512581.190512078.5961808012512581.190612075.451808012512581.190712090.4871808012512581.190812083.8481808012512581.190912073231.398012512581.1901012073198.488012512581.1901112073212.648012512581.1901212073190.558012512581.190131207318075.85712512581.190141207318065.95812512

27、581.190151207318087.25812512581.190161207318097.82412512581.190171207318080150.7112581.190181207318080141.5812581.190191207318080132.3712581.190201207318080156.9312581.190211207318080125138.8881.190221207318080125131.2181.190231207318080125141.7181.190241207318080125149.2981.19025120731808012512560.

28、5829026120731808012512570.9629027120731808012512564.8549028120731808012512575.5299029120731808012512581.1104.8430120731808012512581.1111.2231120731808012512581.198.09232120731808012512581.1120.44表2 各线路的潮流值（各方案与表1相对应，单位：MW）方案线路1234560164.78140.87-144.25119.09135.44157.691165.81140.13-145.14118.63135.

29、37160.762165.51140.25-144.92118.7135.33159.983167.93138.71-146.91117.72135.41166.814166.79139.45-145.92118.13135.41163.645164.94141.5-143.84118.43136.72157.226164.8141.13-144.07118.82136.02157.57165.59143.03-143.16117.24139.66156.598165.21142.28-143.49117.96137.98156.969167.43140.82-152.26129.58132.

30、04153.610165.71140.82-147.08122.85134.21156.2311166.45140.82-149.33125.75133.28155.0912165.23140.85-145.82121.16134.75156.7713164.23140.73-144.18119.12135.57157.214163.04140.34-144.03119.31135.97156.3115165.54141.1-144.32118.84135.06158.2616166.88141.4-144.34118.67134.67159.2817164.07143.03-140.9711

31、8.75133.75158.8318164.27142.29-142.15118.85134.27158.3719164.57141.44-143.3119134.88158.0120163.89143.61-140.25118.64133.28159.1221166.35139.29-144.2119.1136.33157.5922165.54140.14-144.19119.09135.81157.6723166.75138.95-144.17119.15136.55157.5924167.69138.07-144.14119.19137.11157.6525162.21141.21-14

32、4.13116.03135.5154.2626163.54141-144.16117.56135.44155.9327162.7141.14-144.21116.74135.4154.8828164.06140.94-144.18118.24135.4156.6829164.66142.27-147.2120.21135.28157.6530164.7142.94-148.45120.68135.16157.6331164.67141.56-145.88119.68135.29157.6132164.69143.84-150.34121.34135.12157.64表3 各机组的段容量 （单位

33、：MW）机组段1234567891017005000300004023002081562008311004003002040040455510101010150015755150151501010106950102001510200107501551510105103287002002002010155表4 各机组的段价（单位：元/兆瓦小时，记作元/MWh）机组段123456789101-50501241682102523123303634892-56001822032453003203604104953-61001521892332583083564155004-50015017020025

34、53023253804358005-59001161461882152503103965106-60701591732052523053804055207-5001201802512603063153353485488-800153183233253283303318400800表5 各机组的爬坡速率 （单位：MW/分钟）机组12345678速率1.821.41.8表6 各线路的潮流限值（单位：MW）和相对安全裕度线路123456限值165150160155132162安全裕度13%18%9%11%15%14%附录2function k=zhen(fa8,r6) %求k初值

35、 p0=164.78 140.87 -144.25 119.09 135.44 157.69;pj=120 73 180 80 125 125 81.1 90;j=0;for i=1:8 j=fa8(i)-pj(i)+j;endpi=r6-p0;k=pi./j;附录3function ok=zhen1(fa8,r6) %求k系数clcfor i=1:32 a=fa8(i,:);b=r6(i,:); k(i,:)=zhen(a,b);endoptions=;for i=1:8 for j=1:6 ok(i,j)=mean(k(4*i,j),k(4*i-1,j),k(4*i-2,j),k(4*i-

36、3,j); end end附录4function c=zhen4(ok,fa8,r6) %求c值for k=1:32 if k=11 c(k,:)=c(k-1,:); end for j=1:6 c(k,j)=r6(k,j)-fa8(k,:)*ok(:,j); endendc=mean(c);附录5function pout=zhen5(pin) %求潮流值ok=0.08284 0.04828 0.05297 0.11993 -0.02544 0.12201 0.12158 -0.00123;. -0.05456 0.12785 -0.00003 0.03328 -0.08685 -0.112

37、44 -0.01893 -0.09873;. -0.06954 0.06105 -0.15662 -0.00992 0.12449 0.00212 0.00251 -0.20139;. -0.03446 -0.10241 0.20516 -0.02083 -0.01183 -0.00595 0.14492 -0.07655;. 0.00053 0.24329 -0.06455 -0.04113 -0.06522 0.07034 -0.00426 -0.00891;. 0.23781 -0.06017 -0.07787 0.09298 0.04690 0.00008 0.16593 0.0006

38、9;c=110.2965 131.2289 -108.873 77.4817 132.97447 120.6633;pout=ok*pin'+c'附录6function f=zhenfun2(x)x0=120 73 180 80 125 125 81.1 90;x1=150 81 218.2 99.5 135 150 102.1 117;v=2.2 1 3.2 1.3 1.8 2 1.4 1.8;duan=70, 0 ,50 ,0, 0, 30 ,0, 0 ,0, 40;. 30, 0, 20 ,8 ,15 ,6 ,2, 0, 0, 8;. 110 ,0 ,40 ,0, 30

39、,0 ,20, 40, 0, 40;. 55 ,5 ,10, 10, 10, 10, 15 ,0 ,0 ,1;. 75, 5 ,15 ,0 ,15 ,15, 0,10 ,10 ,10;. 95, 0, 10, 20, 0, 15 ,10 ,20 ,0 ,10;. 50 ,15 ,5 ,15 ,10, 10 ,5 ,10 ,3 ,2;. 70, 0, 20, 0 ,20, 0 ,20, 10, 15, 5; price=-505 0 124 168 210 252 312 330 363 489;. -560 0 182 203 245 300 320 360 410 495;. -610 0

40、152 189 233 258 308 356 415 500;. -500 150 170 200 255 302 325 380 435 800;. -590 0 116 146 188 215 250 310 396 510;. -607 0 159 173 205 252 305 380 405 520;. -500 120 180 251 260 306 315 335 348 548;. -800 153 183 233 253 283 303 318 400 800; n0=6 6 5 6 6 6 5 10;b=0 0 0 0 0 0 0 0 ;p=0 0 0 0 0 0 0 0

41、 ;for i=1:8 if x(i)-x0(i)<0 xi=x(i)-x0(i); t=0; if xi<0 for n=n0(i):-1:1 b(i)=duan(i,n)+b(i); if b(i)+xi>=0 b(i)=-duan(i,n)+b(i); break; end end for ni=n0(i):-1:n if ni=n if duan(i,ni)=0 p(i)=-duan(i,ni)2/v(i)/2*(303-price(i,ni)+p(i); t=t+duan(i,ni)/v(i); end else if xi<=x0(i)-x1(i) p(i)

42、=-(xi+b(i)2/v(i)/2*(303-price(i,ni)+(0.25-t-(x1(i)-x0(i)/v(i)*(xi+b(i)*(303-price(i,ni)+p(i); else p(i)=-(xi+b(i)2/v(i)/2*(303-price(i,ni)+(0.25-t-(x1(i)+xi-x0(i)/v(i)*(xi+b(i)*(303-price(i,ni)+p(i); end end end else for n=n0(i):10 b(i)=duan(i,n)+b(i); if b(i)-xi>=0 b(i)=-duan(i,n)+b(i); break; end end for ni=n0(i):n if ni=n if duan(i,ni)=0 p(i)=duan(i,ni)2/v(i)/2*(price(i,ni)-303)+p(i); t=t+duan(i,ni)/v(i); end else p(i)=(xi-b(i)2/v(i)/2*(price(i,ni)-303)+(0.25-t)*(xi-b(i)*(303-price(i,ni)+p(i); end endend elseif x(i)-x1(i)<0 xi=x(i)-x1(i); t=0; if xi<0 for n=n0(i):-1:1 b(