Browse Source

Added code files used in the base case and the computation for the marginal impact of the transmission constraints.

master
Carmen Li 6 months ago
commit
51168a2fec
7 changed files with 5250 additions and 0 deletions
  1. +1
    -0
      VRE/readme.md
  2. +5003
    -0
      output_maximisation/basecase.ipynb
  3. +106
    -0
      output_maximisation/constrained.m
  4. +17
    -0
      output_maximisation/constraint/transmission_constraints.csv
  5. +94
    -0
      output_maximisation/obj_linearP.m
  6. +3
    -0
      output_maximisation/readme.md
  7. +26
    -0
      readme.md

+ 1
- 0
VRE/readme.md View File

@ -0,0 +1 @@
Python code used to download the wind and solar output data.

+ 5003
- 0
output_maximisation/basecase.ipynb
File diff suppressed because it is too large
View File


+ 106
- 0
output_maximisation/constrained.m View File

@ -0,0 +1,106 @@
M_stats=readmatrix('mean.csv');
v_mean = M_stats(:,2);
T_stats=readtable('mean.csv');
zonetechs=T_stats(:,1);
clear M_stats;
clear T_stats;
M_cov = readmatrix('wind_solar_cov_day_night.csv');
M_cov = M_cov(:, [2:end]);
d=32;
% x_i >=0 constraint
lb=zeros(d,1);
ub=ones(d,1);
% equality constraint
AeqI=ones(1,d);
beqI=1;
% fixed power constraint on each INTERIOR points on the frontier
AeqR=v_mean';
Aeq=[ AeqI; AeqR ] ;
% transmission inequality constraints
trans_con=readtable('constraint/transmission_constraints.csv');
trans_con=trans_con(:,2);
trans_con=table2array(trans_con);
M_con=zeros(1,32);
for i =1:d/2
%v_con=cat(2, zeros(1, i-1), power_output(i), zeros(1, 15), power_output(16+i), zeros(1,16-i));
v_con=cat(2, zeros(1, i-1), 1, zeros(1, 15), 1, zeros(1,16-i));
M_con=cat(1, M_con, v_con);
end
M_con=M_con([2:end], :);
% find the lowest risk point
[x_minrisk , fval_minrisk, exit_minrisk, out_minrisk, lambda_minrisk]=quadprog(M_cov,[],M_con,trans_con,AeqI,beqI,lb,ub) ;
S_minrisk= sqrt(2* fval_minrisk) ;
P_minrisk= (v_mean')*x_minrisk ;
% find the highest power point
[x_maxp , mP_maxp, exit_maxp, out_maxp, lambda_maxp]=linprog(-v_mean,M_con,trans_con,AeqI,beqI,lb,ub);
P_maxp=-mP_maxp;
S_maxp= sqrt((x_maxp')*M_cov*x_maxp) ;
% min and max of the sampling points in the interior of the frontier
P_maxp_rd=round(P_maxp , 3);
P_minrisk_rd=round(P_minrisk, 3) ;
if P_minrisk_rd<P_minrisk;
sample_min=P_minrisk_rd+0.0005
else
sample_min= P_minrisk_rd;
end
if P_maxp_rd > P_maxp;
sample_max=P_maxp_rd-0.0005;
else
sample_max=P_maxp_rd;
end
% generate sampling points, more towards the lower end
samples=cat(2,sample_min:0.0005:sample_min+0.0035, sample_min+0.004:0.005: floor(sample_max*100)/100, sample_max, P_maxp-0.000005);
% dummy first entry to store the solutions
v_risk=[0];
M_x=zeros(d,1);
M_lambda_z=zeros(d/2 ,1);
% main loop to find the frontier
for i=samples
beqR = i;
beq=[beqI ; beqR] ;
% find the interior points on the frontier
[x , fval, exit_, out_, lambda_]=quadprog(M_cov,[],M_con,trans_con,Aeq,beq,lb,ub) ;
v_risk=cat(2, v_risk, sqrt(2*fval));
M_x=cat(2, M_x, x);
M_lambda_z=cat(2, M_lambda_z, lambda_.ineqlin);
end
% delete dummy column
v_risk=v_risk(2:end);
M_x=M_x(: , [2:end]);
M_lambda_z=M_lambda_z(: , [2:end]);
% make nice output table
output=cat(2, P_minrisk, samples, P_maxp);
output=cat(1, output, cat(2, S_minrisk, v_risk, S_maxp));
output=cat(1, output, cat(2, x_minrisk, M_x, x_maxp));
lambda_z=cat(2, lambda_minrisk.ineqlin, M_lambda_z, lambda_maxp.ineqlin);
writematrix(lambda_z,'lambda_dR2dXmax.csv') ;
T = array2table(output);
zonetechs=table2cell(zonetechs);
zonetechs=reshape(zonetechs,1,32);
T.Properties.RowNames = [{'mean_P'}, {'mean_s'},zonetechs];
writetable(T,'solution_constrained.csv','WriteRowNames',true);

+ 17
- 0
output_maximisation/constraint/transmission_constraints.csv View File

@ -0,0 +1,17 @@
,X_bar
NSA,0.08163425405989325
ADE,0.1588145857413226
SESA,0.07176232115097092
NQ,0.09257177552582982
CQ,0.11506151427548097
SWQ,0.19362223044297147
SEQ,0.34765823328876133
NNS,0.13186888029783053
NCEN,0.5124289671008763
SWNSW,0.09839742382087051
CAN,0.2526616080396274
CVIC,0.0697678554959203
NVIC,0.13217487570392325
MEL,0.5654584983157979
LV,0.22919710554940376
TAS,0.0674722490472888

+ 94
- 0
output_maximisation/obj_linearP.m View File

@ -0,0 +1,94 @@
M_stats=readmatrix('mean.csv');
v_mean = M_stats(:,2);
T_stats=readtable('mean.csv');
zonetechs=T_stats(:,1);
clear M_stats;
clear T_stats;
M_cov = readmatrix('wind_solar_cov_day_night.csv');
M_cov = M_cov(:, [2:end]);
T_R=readmatrix('solution_constrained.csv');
v_R=T_R(2, (2:end));
v_R(1)=v_R(1)+0.00000001;
v_R(end)=v_R(end)-0.000001;
M_sol_x=T_R((3:end), (2:end));
clear T_R;
d=32;
% x_i >=0 constraint
lb=zeros(d,1);
ub=ones(d,1);
% equality constraint
AeqI=ones(1,d);
beqI=1;
Aeq= AeqI ;
beq=beqI;
% transmission inequality constraints
trans_con=readtable('constraint/transmission_constraints.csv');
trans_con=trans_con(:,2);
trans_con=table2array(trans_con);
M_con=zeros(1,32);
for i =1:d/2
v_con=cat(2, zeros(1, i-1), 1, zeros(1, 15), 1, zeros(1,16-i));
M_con=cat(1, M_con, v_con);
end
M_con=M_con([2:end], :);
%x0=cat(1, zeros(13,1), 1, zeros(18,1));
options = optimset('Algorithm','active-set');
% dummy first entry to store the solutions
v_power=[0];
M_x=zeros(d,1);
M_lambda_z=zeros(d/2 ,1);
j=1;
for i=v_R
fixed_R2= i^2;
x0=M_sol_x(:,j);
c = @(x) [] ;
R2ceq = @(x) x'*M_cov*x - fixed_R2; % fixed R2 quadratic equality constaint
R2eqcon = @(x) deal(c(x),R2ceq(x));
E_p = @(x) -(v_mean'*x);
[x,fval,exitflag,out,lambda,grad,hessian] = fmincon(E_p , x0, M_con, trans_con, Aeq, beq, lb, ub, R2eqcon , options);
v_power=cat(2, v_power, -fval);
M_x=cat(2, M_x, x);
M_lambda_z=cat(2, M_lambda_z, lambda.ineqlin);
j=j+1;
end
% delete dummy column
v_power=v_power(2:end);
M_x=M_x(: , [2:end]);
M_lambda_z=M_lambda_z(: , [2:end]);
% make nice output table
output=cat(1, v_power, v_R, M_x);
T = array2table(output);
zonetechs=table2cell(zonetechs);
zonetechs=reshape(zonetechs,1,32);
T.Properties.RowNames = [{'mean_P'}, {'mean_s'},zonetechs];
writetable(T,'solution_constrained_linobj.csv','WriteRowNames',true);
writematrix(M_lambda_z,'lambda_dPdXmax_linobj.csv') ;

+ 3
- 0
output_maximisation/readme.md View File

@ -0,0 +1,3 @@
The IPython notebook was used to generate the statistical parameters for the output maximising case as well as the analysis plots.
The MATLAB files were used to compute the infinitesimal impact on the linear output and quadratic risk objectives by an infinitesimal increase in the capacity limit in each zone.

+ 26
- 0
readme.md View File

@ -0,0 +1,26 @@
# A Portfolio approach to wind and solar deployment in Australia
## Overview
This repository contains the source codes used to generate and analyse the results of [1]. The directories are named using the same terminology as in the paper.
Wind and solar capacity factor time series are obtained from [2].
Demand data are taken from the "P5MIN_REGIONSOLUTION" dataset published on the AEMO Market Data site [3]. We take the "DEMAND_AND_NONSCHEDGEN" column in the data set as the gross demand and "TOTALDEMAND" as the operational demand. The raw data captured from AEMO are recorded state by state; to disaggregate the state demand into zone demand we use the population proportion method illustrated in [4].
The current transmission
limits between zones are also estimated using the data from Xenophon and Hill (2018).
## References
[1] C. K. Chong, C. Li, D. Reiner and F. Roques, A Portfolio approach to wind and solar deployment in Australia. EPRG working papers 2022 (2020), at [https://www.eprg.group.cam.ac.uk/eprg-working-paper-2022/](https://www.eprg.group.cam.ac.uk/eprg-working-paper-2022/)
[2] S. Pfenninger and I. Staffell Renewables.ninja, at [https://www.renewables.ninja/](https://www.renewables.ninja/)
[3] Australian Energy Market Operator (AEMO). NEMWEB Market Data, at [http://www.nemweb.com.au/] (http://www.nemweb.com.au/)

Loading…
Cancel
Save