Add files via upload

main_Figure6.m

@@ -0,0 +1,146 @@
+clc
+clear all
+close all
+%% This script allows drawing the figure 6
+% Use Section 1 to estimate the energy transfer as a function
+% of ageing limit. (data generation
+
+% Otherwise, you can dircetly launch section 2 to construct the figure "ageing 
+% limit vs Energy transfer" from precalculated data.('ONAN_interm_results.mat')
+
+% Section 3 draws the figure: Energy transfer as a function of time
+
+% Section 4 prepares the data for Figure 8 
+%% Section 1: Generation of data
+% Initial data 
+TIM=linspace(1,1440,1440)'; % time in minutes for 1 day
+Temperature=-50:1:50;       % Full range of ambient temperatures
+Ageing_limit=0.1:0.1:12.7;  % Range of Ageing_limits 
+AMB=linspace(20,20,24)';    % Fixed ambient temperature, hour time step
+
+% Finding the energy transfer for each ageing limit
+for i=1:length(Ageing_limit) % for ageing limit
+
+    % Solve the optimization problem for 1  day
+    [Energy_limit,HST_limit,AEQ_limit,TOT_limit,Energy]=optimal_energy_limitFig6(AMB,Ageing_limit(i));
+
+    % Save energy transfer. So far, units are  pu*1min but later we will
+    % normalize them
+    Energy_result(i)=Energy;
+
+    % Energy limit (loading profile of transformer in pu)
+    PUL_result(1:1440,i)=Energy_limit;
+
+    % Loss of Life results
+    LoL_result(i)=AEQ_limit;
+
+    % caclulate the ageing acceleration factor (DL)
+    for m=1:length(HST_limit)
+        DL(m,:) = (2^((HST_limit(m)-98)/6));
+    end
+
+    % Calculate value of cummulated ageing at each moment
+    Current_ageing=0;
+    for j=1:length(DL)
+        Current_ageing(j)=Current_ageing(end)+DL(j);
+    end
+
+    %Save the result of cummulated ageing for given ageing limit (i)
+    Current_ageing_result(1:1440,i)=Current_ageing;
+    i % show the iteration
+end
+%% Section 2: Construct the figure ageing limit vs Energy transfer 
+clear;clc;close all
+% Load the results from previous section
+load('ONAN_interm_results.mat')
+
+
+% Choose the calendar life left (see cases in Table 1 of the paper)
+studied_case=1; % 1,2 or 3
+
+if studied_case==1
+    % I case
+    calendar_time_left=50;   % 50%
+    insulation_time_left=70; % 70%
+
+elseif studied_case==2
+    % II case
+    calendar_time_left=50;   % 50%
+    insulation_time_left=50; % 70%
+
+elseif studied_case==3
+    % III case
+    calendar_time_left=50;   % 50%
+    insulation_time_left=30; % 70%
+else
+    error('Check the value of studied_case. It shoul be 1,2 or 3')
+end
+
+% Find the ratio (optimal ageing limit)
+ratio=insulation_time_left/calendar_time_left;
+
+% Find the calendar life in minutes of days
+calendar_time_left_min=calendar_time_left*1440/100;
+
+% Find the index which shows how fast insulation life left will be consumed at
+% given ageing limit
+for i=1:length(Ageing_limit)
+    index(i)=(insulation_time_left*1440/100)/Ageing_limit(i);
+end
+
+% Round the index to get an integer (approximation!)
+index=round(index);
+
+% Find time for which insulation resource can be used. 
+time_percent=(index./1440); % in pu
+% Note that we calculate it relatively to 1 day 1440. Units are pu!
+
+% Find how much energy is transfered for each ageing limit
+for i=1:length(Ageing_limit)
+    if index(i)>calendar_time_left_min
+        ind=calendar_time_left_min;
+    else
+        ind=index(i);
+    end
+    Energy_day(i)=ind/calendar_time_left_min*sum(PUL_result(:,i));
+end
+
+% Energy_day=Energy_result;
+
+% Normalize the energy transfer relatively to ageing limit =1 pu
+Energy_day_percent=Energy_day/Energy_day(10)*100;
+
+plot(Ageing_limit,Energy_day_percent)
+ylabel('Energy, % of Energy at Ageing limit=1')
+xlabel('Ageing limit(AAF), pu')
+
+% time_percent=index/1440*100;
+%% Section 3:Energy transfer as a function of time (additonal figure)
+% Assumption: we assume that for each day the loading profile (and thus energy
+% transfer) are the same. Thus, we can estimate how much energy can be
+% transfered during the entire life if we multiply the energy trasnfer for
+% 1 day on possible operation at given ageing limit (in per units!)
+Energy_all_life=time_percent.*Energy_day;
+
+% Normalize the energy transfer for each ageing limit (relatively to ageing
+% limit = 1pu)
+Energy_all_life=Energy_all_life./Energy_all_life(10)*100;
+
+% Convert time_percent from pu into percent
+time_percent=time_percent*100;
+
+% Set the intitial position 
+start=[0 0];
+
+figure % create the figure
+hold on 
+
+% Draw a figure: Energy transfer as a function of time (for each ageing
+% limit)
+for i=1:length(Ageing_limit)
+    name=['Ageing limit ',num2str(Ageing_limit(i))];
+   plot([50 50+time_percent(i)], [30 30+Energy_all_life(i)],'DisplayName',name)
+end
+
+ylabel('Energy, % of Energy at Ageing limit=1')
+xlabel('Time, % of Time at Ageing limit=1')

optimal_energy_limitFig6.m

@@ -0,0 +1,48 @@
+function [Energy_limit,HST_limit,AEQ_limit,TOT_limit,Energy,exitflag]=optimal_energy_limit(AMB,AAF_limit)
+%% Problem-based for DTR transformer
+nHours = numel(AMB);
+%% Define the optimization problem and the optimization variables
+DTRprob=optimproblem('ObjectiveSense','maximize'); 
+% DTR in an hour for transformer
+DTR = optimvar('DTR',nHours,'LowerBound',0,'UpperBound',1.5);
+%% Define the objective function
+% Energy_transfer maximization;
+Energy_transfer = sum(DTR);
+
+% showexpr(Energy_transfer)
+% set objective
+DTRprob.Objective = Energy_transfer;
+
+%% DTR constraints
+%type OF
+[HST_max,TOT_max,~,~,AEQ,~]=fcn2optimexpr(@distrbution_transformer_optim,DTR,AMB);
+% Ageing below or equal normal
+DTRprob.Constraints.ageing=AEQ<=AAF_limit;
+% HST below or equal 120 
+DTRprob.Constraints.hst=HST_max<=120;
+% TOT below or equal to 105
+DTRprob.Constraints.tot=TOT_max<=105;
+% HST start equal to HST end
+% DTRprob.Constraints.hst_equality1=HST_1==HST_end;
+% showconstr(DTRprob.Constraints.ageing)
+
+% showproblem(DTRprob)
+%% OPtimal solution
+
+% First apprximation
+% [DTR_approx]=method1_IEC_improved(AMB);
+DTR_approx=linspace(1,1,length(AMB))';
+x0.DTR=DTR_approx';
+
+% call the optimization solver to find the best solution
+
+% [sol,TotalCost,exitflag,output] = solve(DTRprob,x0,'Options',options);
+options=optimset('disp','iter','LargeScale','off','TolFun',.001,'MaxIter',100000000,'MaxFunEvals',1000000000);
+options.Algorithm='sqp';
+[sol,~,exitflag,~] = solve(DTRprob,x0,'Options',options);
+%% Energy comparison;
+Energy_limit=sol.DTR;
+Energy=sum(sol.DTR);
+Energy_limit=PUL_to_1min(Energy_limit,60);
+[HST_limit,TOT_limit,AEQ_limit,~,~]=distrbution_transformer_full(Energy_limit,AMB);
+end