env.out()<<\ cout<<\ }
catch(IloException& e){ cerr<<\ }
catch(...){ cerr<<\ }
env.end(); return 0; }
运行结果如下图所示：
同样的，上述cplex程序仅进行了对3个小时内最优解的求解。对半天12小时我们依然采用了较为简单的matlab编程,所编程序如下：
P1=sdpvar (1,12); %P1代表1号机组12小时内每小时的功率 P2=sdpvar (1,12); %P2代表2号机组12小时内每小时的功率 P3=sdpvar (1,12); %P3代表3号机组12小时内每小时的功率 %定义每台机组每小时功率P A=P1; B=P2;
C=P3;%将功率赋予数组
s1=binvar (1,12);%S1代表1号机组12小时内每小时的开关状态 s2=binvar (1,12);%S2代表2号机组12小时内每小时的开关状态 s3=binvar (1,12);%S3代表3号机组12小时内每小时的开关状态 %定义0-1变量 S1=s1; S2=s2;
11
S3=s3;%将0-1变量赋予数组
S=[1 1 1 1 1 1 1 1 1 1 1 1];
Pw1=30;%Pw1代表机组1的爬坡功率最大值 Pw2=50;%Pw2代表机组2的爬坡功率最大值 Pw3=70;%Pw3代表机组3的爬坡功率最大值 %定义每台机组每小时功率最大改变量Pw
f=126+0.3*A*S'+0.2*B*S'+0.25*C*S';%求最小开销
F=set(10<=A)+set(A<=100)+set(20<=B)+set(B<=368)+set(30<=C)+set(C<=389);%P约束
F=F+set(A(1)+B(1)+C(1)==87)+set(A(2)+B(2)+C(2)==157)+set(A(3)+B(3)+C(3)==267)+set(A(4)+B(4)+C(4)==355);
F=F+set(A(5)+B(5)+C(5)==399)+set(A(6)+B(6)+C(6)==402)+set(A(7)+B(7)+C(7)==379)+set(A(8)+B(8)+C(8)==346);
F=F+set(A(9)+B(9)+C(9)==318)+set(A(10)+B(10)+C(10)==302)+set(A(11)+B(11)+C(11)==303)+set(A(12)+B(12)+C(12)==318); %任意小时功率P之和约束
F=F+set(abs(A(2)-A(1))<=Pw1); F=F+set(abs(A(3)-A(2))<=Pw1); F=F+set(abs(A(4)-A(3))<=Pw1); F=F+set(abs(A(5)-A(4))<=Pw1); F=F+set(abs(A(6)-A(5))<=Pw1); F=F+set(abs(A(7)-A(6))<=Pw1); F=F+set(abs(A(8)-A(7))<=Pw1); F=F+set(abs(A(9)-A(8))<=Pw1); F=F+set(abs(A(10)-A(9))<=Pw1); F=F+set(abs(A(11)-A(10))<=Pw1); F=F+set(abs(A(12)-A(11))<=Pw1); %Pw1爬坡约束
F=F+set(abs(B(2)-B(1))<=Pw2); F=F+set(abs(B(3)-B(2))<=Pw2); F=F+set(abs(B(4)-B(3))<=Pw2); F=F+set(abs(B(5)-B(4))<=Pw2); F=F+set(abs(B(6)-B(5))<=Pw2); F=F+set(abs(B(7)-B(6))<=Pw2); F=F+set(abs(B(8)-B(7))<=Pw2); F=F+set(abs(B(9)-B(8))<=Pw2); F=F+set(abs(B(10)-B(9))<=Pw2); F=F+set(abs(B(11)-B(10))<=Pw2); F=F+set(abs(B(12)-B(11))<=Pw2); %Pw2爬坡约束
F=F+set(abs(C(2)-C(1))<=Pw3); F=F+set(abs(C(3)-C(2))<=Pw3); F=F+set(abs(C(4)-C(3))<=Pw3); F=F+set(abs(C(5)-C(4))<=Pw3);
12
F=F+set(abs(C(6)-C(5))<=Pw3); F=F+set(abs(C(7)-C(6))<=Pw3); F=F+set(abs(C(8)-C(7))<=Pw3); F=F+set(abs(C(9)-C(8))<=Pw3); F=F+set(abs(C(10)-C(9))<=Pw3); F=F+set(abs(C(11)-C(10))<=Pw3); F=F+set(abs(C(12)-C(11))<=Pw3); %Pw3爬坡约束
solvesdp(F,f);%进行非线性规划 D=[A,B,C];
E=reshape(D,12,3)'; double(f);
double(E);%输出结果
运行结果如下：
P1=[10 10.0000000000000 10 10 10.0000000000000 10.0000000000000 9.99999999999999 10 10 10.0000000000000 10 10.0000000000000]
P2=[47.0000000000000 97.0000000000000 147.000000000000 197.000000000000 247.000000000000 297.000000000000 339.000000000000 306.000000000000 278.000000000000 262.000000000000 263.000000000000 278.000000000000] P3=[30.0000000000000 49.9999999999999 110.000000000000 148.000000000000 142.000000000000 94.9999999999998 30 30 30 30 30 30] Min f=902.350000000000
以上我们仅对12小时内的情况进行求解，因为在执行24小时的优化过程中，会出现因为约束过多或数值不符而导致没有可行域的结果。因此我们采用了12小时的优化方法，对12-24小时采用动态规划的方法，得出最优解，方法同上。
３.４问题3案例分析：
对于问题3，模型已变为非线性规划，而我们对cplex还不甚熟悉，不知道如何运用cplex解决非线性规划问题，所以我们只采用了matlab进行编程。
由于我们模型的决策变量之于24小时高达144个，起初运行时不仅耗费时间，而且得到的最优解无法应用于实际情况，因此我们采用6小时决策，即每6小时进行一次规划，得到最优解。例如，对1-6小时编写程序如下：
p1=sdpvar (1,6); p2=sdpvar (1,6);
p3=sdpvar (1,6);%定义每台机组每小时功率p PA=sdpvar (1,6); PB=sdpvar (1,6); PC=sdpvar (1,6); A=p1; B=p2;
C=p3;%将功率赋予数组 s1=binvar (1,6); s2=binvar (1,6);
s3=binvar (1,6);%定义0-1变量
13
S1=s1; S2=s2;
S3=s3;%将0-1变量赋予数组 PA(1)=A(1)*S1(1); PA(2)=A(2)*S1(2); PA(3)=A(3)*S1(3); PA(4)=A(4)*S1(4); PA(5)=A(5)*S1(5); PA(6)=A(6)*S1(6); PB(1)=B(1)*S2(1); PB(2)=B(2)*S2(2); PB(3)=B(3)*S2(3); PB(4)=B(4)*S2(4); PB(5)=B(5)*S2(5); PB(6)=B(6)*S2(6); PC(1)=C(1)*S3(1); PC(2)=C(2)*S3(2); PC(3)=C(3)*S3(3); PC(4)=C(4)*S3(4); PC(5)=C(5)*S3(5); PC(6)=C(6)*S3(6); pw1=30; pw2=50;
pw3=70;%定义每台机组每小时功率最大改变量pw
f=126+0.3*A*S1'+0.2*B*S2'+0.25*C*S3';%求最小开销
F=set(20<=A)+set(A<=100)+set(50<=B)+set(B<=368)+set(75<=C)+set(C<=389);%p约束
F=F+set(PA(1)+PB(1)+PC(1)==87)+set(PA(2)+PB(2)+PC(2)==157)+set(PA(3)+PB(3)+PC(3)==267)+set(PA(4)+PB(4)+PC(4)==355);
F=F+set(PA(5)+PB(5)+PC(5)==399)+set(PA(6)+PB(6)+PC(6)==402);%任意小时功率p之和约束
F=F+set(S1(1)*S1(2)*abs(PA(1)-PA(2))<=pw1); F=F+set(S1(2)*S1(3)*abs(PA(2)-PA(3))<=pw1); F=F+set(S1(3)*S1(4)*abs(PA(3)-PA(4))<=pw1); F=F+set(S1(4)*S1(5)*abs(PA(4)-PA(5))<=pw1); F=F+set(S1(5)*S1(6)*abs(PA(5)-PA(6))<=pw1); %pw1爬坡约束
F=F+set(S2(1)*S2(2)*abs(PB(1)-PB(2))<=pw2); F=F+set(S2(2)*S2(3)*abs(PB(2)-PB(3))<=pw2); F=F+set(S2(3)*S2(4)*abs(PB(3)-PB(4))<=pw2); F=F+set(S2(4)*S2(5)*abs(PB(4)-PB(5))<=pw2); F=F+set(S2(5)*S2(6)*abs(PB(5)-PB(6))<=pw2); %pw2爬坡约束
F=F+set(S3(1)*S3(2)*abs(PC(1)-PC(2))<=pw3);
14
F=F+set(S3(2)*S3(3)*abs(PC(2)-PC(3))<=pw3); F=F+set(S3(3)*S3(4)*abs(PC(3)-PC(4))<=pw3); F=F+set(S3(4)*S3(5)*abs(PC(4)-PC(5))<=pw3); F=F+set(S3(5)*S3(6)*abs(PC(5)-PC(6))<=pw3); %pw3爬坡约束
solvesdp(F,f);%进行非线性规划 D=[PA,PB,PC,S1,S2,S3]; E=reshape(D,6,6)'; double(f);
double(E);%输出结果
运行结果如下：
P1 34.32028145 59.64081816 89.86014999 38 32 P2 165.6426186 258.6932121 267 317 367 P3 87 157 206.128554 276.3176527 388.130609 S1 0 0 0 1 1 S2 0 0 1 1 1 S3 1 1 0 0 0
Min f=482
相同的，在6-12,12-18,18-24小时可以在1-6小时基础上进行动态规划，从而得到在另外三个时段的最优解，方法同问题3。
四、研究结论：
由三个问题的建模和编程我们发现，当一个问题是比较简单的线性规划问题的时候，利用软件去解决是比较方便的，即使变量数目很多，在很短的时间内也可以解决（如果有最优解的话），但是当问题变得复杂时，不仅仅变量和约束变多，问题也很有可能转化为非线性规划问题，这样再利用软件求解时，不仅程序运行的时间长，而且很可能得不到最优解或者得到的最优解不能运用于实际情况。在这种情况下，我们对此有两种方案：1.重新对问题进行建模，对模型进行优化，使之能够尽量的用线性约束表示出来，减小软件实现的难度和所需要的时间；2.如果模型已足够优化，但问题的模型还是多变量非线性规划，我们通过减少决策变量的方法来先进行小部分决策变量的最优解的求解，在此基础上进行另外一小部分决策变量最优解的求解，依次类推，最终实现全变量的最优解的求解。这种方法实际上是一种动态规划的方法，是在我们对全部决策变量寻找最优解不可行的情况下得到的最优解，并不是真正的“最优解”，而是一种次优解，但这对我们的现实情况是具有指导意义的。
对于这个数学模型，如果我们从更深的一方面去思考，就会发现更多的东西。首先，这种对机组的调度仅仅建立在发电成本的基础上，而这在现实中肯定是不够合理的。在现实中，当我们考虑一个发电机组的动态调度问题时，环境污染，启停时间限制，电能的运输损耗，人工成本，以及发电边际成本的最优取值等都是不可忽略的因素。如果在这个建模问题中，老师没有给出一个明确标准，还是让学生们自己在考虑不同因素的基础上给出一个“最好的”的动态调度方案，那么，每个组给出的评价标
15