Nunits = 3;
Horizon = 48;
Pmax = [100;50;25];
Pmin = [20;40;1];
Q = diag([.04 .01 .02]);
C = [10 20 20];
Pforecast = 100 + 50*sin((1:Horizon)*2*pi/24);
onoff = binvar(Nunits,Horizon,'full');
P = sdpvar(Nunits,Horizon,'full');
Constraints = [];
for k = 1:Horizon
Constraints = [Constraints, onoff(:,k).*Pmin<= P(:,k) <= onoff(:,k).*Pmax];
end
for k = 1:Horizon
Constraints = [Constraints, sum(P(:,k)) >= Pforecast(k)];
end
Objective = 0;
for k = 1:Horizon
Objective = Objective + P(:,k)'*Q*P(:,k) + C*P(:,k);
end
ops = sdpsettings('verbose',1,'debug',1);
optimize(Constraints,Objective,ops)
stairs(value(P)');
legend('Unit 1','Unit 2','Unit 3');
minup = [6;30;1];
mindown = [3;6;3];
for k = 2:Horizon
for unit = 1:Nunits
% indicator will be 1 only when switched on
indicator = onoff(unit,k)-onoff(unit,k-1);
range = k:min(Horizon,k+minup(unit)-1);
% Constraints will be redundant unless indicator = 1
Constraints = [Constraints, onoff(unit,range) >= indicator];
end
end
for k = 2:Horizon
for unit = 1:Nunits
% indicator will be 1 only when switched off
indicator = onoff(unit,k-1)-onoff(unit,k);
range = k:min(Horizon,k+mindown(unit)-1);
% Constraints will be redundant unless indicator = 1
Constraints = [Constraints, onoff(unit,range) <= 1-indicator];
end
end
ops = sdpsettings('verbose',2,'debug',1);
optimize(Constraints,Objective,ops);
stairs(value(P)');
legend('Unit 1','Unit 2','Unit 3');
Unit3Levels = [0 1 6 10 12 20];
for k = 1:Horizon
Constraints = [Constraints, ismember(P(3,k),Unit3Levels)];
end
optimize(Constraints,Objective);
stairs(value(P)');
legend('Unit 1','Unit 2','Unit 3');
Nhist = max([minup;mindown]);
Pforecast = sdpvar(1,Horizon);
HistoryOnOff = sdpvar(Nunits,Nhist,'full');
DemandSlack = sdpvar(1,Horizon);
PreviusP = sdpvar(3,1);
DemandPenalty = 1000;
ChangePenalty = 1000;
Constraints = [];
Objective = 0;
for k = 1:Horizon
Constraints = [Constraints, onoff(:,k).*Pmin <= P(:,k) <= onoff(:,k).*Pmax];
Constraints = [Constraints, sum(P(:,k))+DemandSlack(k) >= Pforecast(k)];
Constraints = [Constraints, DemandSlack(k) >= 0];
Objective = Objective + P(:,k)'*Q*P(:,k) + C*P(:,k);
Objective = Objective + DemandPenalty*DemandSlack(k);
end
function C = consequtiveON(x,minup)
if min(size(x))==1
x = x(:)';
end
if size(x,1) ~= size(minup,1)
error('MINUP should have as many rows as X');
end
Horizon = size(x,2);
C = [];
for k = 2:size(x,2)
for unit = 1:size(x,1)
% indicator will be 1 only when switched on
indicator = x(unit,k)-x(unit,k-1);
range = k:min(Horizon,k+minup(unit)-1);
% Constraints will be redundant unless indicator = 1
affected = x(unit,range);
if strcmp(class(affected),'sdpvar')
% ISA behaves buggy, hence we use class+strcmp
C = [C, affected >= indicator];
end
end
end
Constraints = [Constraints, mindup,consequtiveON(1-[HistoryOnOff onoff])];
Constraints = [Constraints, mindown,consequtiveON(1-[HistoryOnOff onoff])];
for k = 2:Horizon
Objective = Objective + ChangePenalty*norm(P(:,k)-P(:,k-1),1);
end
Objective = Objective + ChangePenalty*norm(P(:,1)-PreviusP,1);
Parameters = {HistoryOnOff, Pforecast, PreviusP};
Outputs = {P,onoff};
ops = sdpsettings('verbose',2,'debug',1);
Controller = optimizer(Constraints,Objective,ops,Parameters,Outputs);
oldOnOff = repmat([1;1;0],1,Nhist);
oldP = repmat([100;40;0],1,Nhist);
for k = 1:500
% Base-line forecast
forecast = 100;
% Daily fluctuation
forecast = forecast + 50*sin((k:k+Horizon-1)*2*pi/24);
% Some other effect
forecast = forecast + 20*sin((k:k+Horizon-1)*2*pi/24/7);
% and yet some other effect
forecast = forecast + randn(1,Horizon)*5;
[solution,problem] = Controller{oldOnOff, forecast, oldP(:,end)};
optimalP = solution{1};
optimalOnOff = solution{2};
hold off
stairs([-Nhist+1:0 1:Horizon],[oldP optimalP]');
hold on
stairs(1:Horizon,forecast,'k+');
axis([-Nhist Horizon 0 170]);
drawnow
pause(0.05);
% Shift history
oldP = [oldP(:,2:end) optimalP(:,1)];
oldOnOff = [oldOnOff(:,2:end) optimalOnOff(:,1)];
end
end
上述代码有没有错误，如果有请帮我修改
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
144
