unitcommitment.m

%% Basic model
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)
clf
stairs(value(P)');
legend('Unit 1','Unit 2','Unit 3');



%% Adding minimum up- and down-time
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');

%% Quantized power-levels
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');


%% Efficient simulation
Nhist = max([minup;mindown]);
Pforecast = sdpvar(1,Horizon);
HistoryOnOff = sdpvar(Nunits,Nhist,'full');
DemandSlack = sdpvar(1,Horizon);
PreviusP = sdpvar(3,1);
DemandPenalty = 1000;
ChangePenalty = 100;
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
Constraints = [Constraints, consequtiveON([HistoryOnOff onoff],minup)];
Constraints = [Constraints, consequtiveON(1-[HistoryOnOff onoff],mindown)];
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
k
    % 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

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
end