function Results = Master(History,UB,tau)
%% Initialization

[~,k] = size(History);
load('Datos.mat');

for n = 1:Datos.N
    % Coefficients to linearize conic constraints
    phi = deg2rad(360);
    sp(n) = sin((phi*n)/Datos.N) - sin((phi*(n-1))/Datos.N);
    sq(n) = cos((phi*n)/Datos.N) - cos((phi*(n-1))/Datos.N);
    sF(n) = sin(phi/Datos.N);
end

%% Problem creation

DN = optimproblem('ObjectiveSense','minimize');

%% Variables

% Community
pi = optimvar('pi',1,Datos.nt,'LowerBound',0,'UpperBound',Datos.Smax); % Active power imported (kW)
pe = optimvar('pe',1,Datos.nt,'LowerBound',0,'UpperBound',Datos.Smax); % Active power exported (kW)
qi = optimvar('qi',1,Datos.nt,'LowerBound',0,'UpperBound',Datos.Smax); % Reactive power imported (kvar)
qe = optimvar('qe',1,Datos.nt,'LowerBound',0,'UpperBound',Datos.Smax); % Reactive power exported (kvar)
fp = optimvar('fp',Datos.nb,Datos.nt); % Active power flows (kW)
fq = optimvar('fq',Datos.nb,Datos.nt); % Reactive power flows (kvar)
V = optimvar('V',Datos.nb,Datos.nt,'LowerBound',Datos.Vmin,'UpperBound',Datos.Vmax); % Nodal voltage (V)

% Prosumers
pC = optimvar('pC',Datos.npc,Datos.nt,'LowerBound',-Datos.pP,'UpperBound',Datos.pP); % Prosumers net power (kW)

% Auxiliary variable
beta = optimvar('beta',1,Datos.npc,'LowerBound',0);

%% Objective function

Bundle = 0;
for c = 1:Datos.npc
    Bundle = Bundle + (tau(c)/2).*sum(sum((pC(c,:) - History(c,k).pC).^2));
end

DN.Objective = pi*Datos.TOU' - sum(Datos.cex.*pe) + sum(beta) + Bundle;

%% Constraints

% Power flows
DN.Constraints.Activa_SUB = pi - pe == fp(1,:); % Active power flow balance at root node
DN.Constraints.Reactiva_SUB = qi - qe == fq(1,:); % Reactive power flow balance at root node
for b = 1:Datos.nb % Power flow through branches
    FlowsP(b,1:Datos.nt) = sum(pC(Datos.pcb == b,:),1) + sum(fp(Datos.UW == b,:),1) == fp(b,:);
    FlowsQ(b,1:Datos.nt) = sum(pC(Datos.pcb == b,:),1).*Datos.phir + sum(fq(Datos.UW == b,:),1) == fq(b,:);
end
DN.Constraints.FlowsP = FlowsP;
DN.Constraints.FlowsQ = FlowsQ;

% Voltages
Nodal_Voltage(1,1:Datos.nt) = V(1,:) == Datos.V0;
for b = 2:Datos.nb
    Nodal_Voltage(b,1:Datos.nt) = V(b,:) == V(Datos.UW(b),:) - 1e3.*(fp(b,:).*Datos.Re(b-1) - fq(b,:).*Datos.X(b-1))/Datos.V0;
end
DN.Constraints.Nodal_Voltage = Nodal_Voltage;

% Power flow limits (linearized)
for n = 1:Datos.N
    for b = 1:Datos.nb
        Flow_Limit(b,n,1:Datos.nt) = sp(n).*fp(b,:) - sq(n).*fq(b,:) - sF(n).*Datos.Smax <= 0;
    end
end
DN.Constraints.Flow_Limit = Flow_Limit;

% Limits on costs of prosumers
% DN.Constraints.UB = sum(beta) <= UB;

%% Cuts

for j = 1:k
    for c = 1:Datos.npc
        C(c,j) = beta(c) >= History(c,j).f - History(c,j).Sen*(pC(c,:) - History(c,j).pC)';
    end
end
          
DN.Constraints.BETA = C;

%% Solver calling

[Results.sol, Results.f,~,~,Results.Sen] = solve(DN);

%% Results of interest

Results.Community_cost = Results.sol.pi*Datos.TOU' - sum(Datos.cex.*Results.sol.pe); % Total community cost (€)
Results.Total_Imported = sum(Results.sol.pi); % Total energy imported by the community (kWh)
Results.Total_Exported = sum(Results.sol.pe); % Total energy exported by the community (kWh)

end