load('input_data.mat')

%% Create the optimization problem 

energyprob = optimproblem;

solar = optimvar('solar',2,'Type','integer','LowerBound',0); % integer because we only want integer values of panels
wind = optimvar('wind',2,'Type','integer','LowerBound',0);

costPanel = [96000, 265000]; % cost for each solar system, $
powerPanel = [30000, 100000]; % power output for each array, Watts

costTurbine = [32899, 85000]; % cost for each turbine, $ 
powerTurbine = [5000, 20000]; % power output for each turbine, Watts

costSolar = costPanel*solar;    % cost times number of arrays
costWind = costTurbine*wind;    % cost times number of turbines
cost = costSolar + costWind;    % total cost of new additions
energyprob.Objective = cost;    % make the objective function the cost

% baseline before additions
peakWind = 36e3;
peakPV = 63e3;  
peakCogen = 1.8e6;

%% Can peak power produced on worst case day meet peak demand? 

peakAdd = (powerPanel*solar + powerTurbine*wind) >= (3e6 - (peakPV + peakCogen + peakWind)); % peak that added solar and wind can produce 
cost = costSolar + costWind <= 12000000;

%% Constrain and solve

energyprob.Constraints.conpl = peakAdd; % the additions plus what is there now must be able to produce 2.5MW peak load 
energyprob.Constraints.conc = cost; % cost less than $12 mil

[sol,fval] = solve(energyprob);
sol.solar
sol.wind 

%% Scale Solar & Wind Power Curves to Include Additions

% start with 1 5kW turbine as the output of the above optimization
% find that it can't meet load all day, scale up using A and B

A = 40;
B = 79;

addSolar = sol.solar(1)*powerPanel(1) + sol.solar(2)*powerPanel(2);        % calculate added solar power
addWind = sol.wind(1)*A*powerTurbine(1) + B*powerTurbine(2);               % calculate added wind power, don't multiply sol.wind(2) because it is zero

% ((total power being added / current peak of power curve) * power curve) + power curve 
% this creates the same shape curve with peak of added power, then sum with what we already have

newSolarPower = ((addSolar/(max(solarPower))))*solarPower + solarPower;     % on cloudy day
newWindPower = ((addWind)/(max(windPower)))*windPower + windPower;          % on slow wind speed day

cloudy = 0.8;   % scale down by 20 percent on worst case day
slow = 0.8;     % scale down by 20 percent on worst case day
cloudySolarPower = (cloudy*((addSolar/(max(solarPower))))*solarPower) + solarPower; % on cloudy day
slowWindPower = (slow*((addWind)/(max(windPower)))*windPower) + windPower;          % on slow wind speed day

totalSupply = cloudySolarPower + slowWindPower + cogenPower; % power generated after additions
solarCost = sol.solar(1)*costPanel(1) + sol.solar(2)*costPanel(2);
windCost = A*sol.wind(1)*costTurbine(1) + B*sol.wind(2)*costTurbine(2);
totalCost = solarCost + windCost

if totalSupply >= Demand
    fprintf('GOOD')
else
    fprintf('Need more power')
end

figure(1)
plot(totalSupply,'b','LineWidth',1.7)
hold on
plot(Demand,'r','LineWidth',1.7)
title('Adding Solar and Wind','fontSize',18)
xlabel('time (seconds)','fontSize',14)
ylabel('MW','fontSize',14)
xlim([0 1441])
legend('Supply','Demand','fontSize',14)
 