Edit on GitHub

Variable DC-link voltages

Since there are not only ideal sources in the regular grid, alternative, variable DC-link voltages shall be presented here. The models to be used will be discussed in detail and it will be explained how they can be created and used in the Env.

  • Modelling a PV-Array,

  • Modelling a Battery,

  • Access the models.

using ElectricGrid
using PlotlyJS
[ Info: Precompiling ElectricGrid [d6749236-802c-4af3-9b09-2fcebf7269d7]

Modelling a PV-Array

This notebook will explain t arrays structurally and how they were implemented.

As a basis serves the paper which explains the physical characteristics of a PV module in more detail. In order to avoid further calculations at runtime, a simpler approximation is used, the ideal single diode model. This neglects the parallel and series resistance of the PV module.

Equivalent curcuit diagram of the PV array

In the following we define a mutable struct with the most important data of a PV module, most of the data are given by the data sheets.

ModulePV = solar_module(); # For this example we use the default values

As a rule, many of the modules are initially connected in series, especially in large systems, in order to increase the voltage. Further modules can then be added in parallel to further increase the current.

We therefore define another mutable struc, which inherits the properties of the PV module and specifies how many modules are to be created in series and then again in parallel to it.

For real applications, this often leads to problems, since shadowing of individual cells (!!!) can result in a large drop in power. Therefore, in practice, countermeasures are often taken, such as the integration of bypass diodes.

ArrayPV = solar_array(;solar_module=ModulePV, serial=1, parallel=1)
solar_array(solar_module(2.0381e-10, 1.2, 1.3806e-23, 1.6022e-19, 273.15, 36, 1000, 0.0013, 3.11), 1, 1, 0.0, 0.0)

Now let us take a look characteristics for different irradiations. Therefore, we define a function that gives us the current as a function of voltage, irradiation and temperature. In addition, a function should also output the voltage as a function of current, irradiation and temperature, since this will later be interesting for our application in ElectricGrid.jl.

Function for module

For the purpose of this notebook, the functions here are defined vectorially so that they can be evaluated directly for multiple irradiances and voltages.

The functions that are later used in the ElectricGrid.jl environment return only a scalar value.

function GetI_vec(SolarArr::solar_array, V, G, T)
    self = SolarArr.solar_module
    
    function I_photo(self::solar_module, G, T)
        dT = self.T_0 + T
        I_ph = G./self.G_ref*(self.I_ph_ref + self.mu_sc * dT).* ones(length(V))'
        return I_ph
    end;

    function I_diode(self::solar_module, V, G, T)
        dT = self.T_0 + T
        V_T = self.k*dT/self.q
        I_d = self.I_0 * (exp.(V  / (self.ni * self.N_cell * SolarArr.serial * V_T)).-1).*ones(5)'
        return I_d
    end;
    
    I = (I_photo(self, G, T)' - I_diode(self, V, G, T)) * SolarArr.parallel
    return I
end

GetI_vec (generic function with 1 method)

Characteristics for PV modules

Here the values are defined for which the modules are to be evaluated.

T = 25
G = collect(200:200:1000)
V = collect(0:0.1:30);
i = GetI_vec(ArrayPV, V, G, T);
labels = ["200 W/m^2" "400 W/m^2" "600 W/m^2" "800 W/m^2" "1000 W/m^2"];
plot([scatter(x=V, y=i[:,1], mode="lines", name=labels[1]),
scatter(x=V, y=i[:,2], mode="lines", name=labels[2]),
scatter(x=V, y=i[:,3], mode="lines", name=labels[3]),
scatter(x=V, y=i[:,4], mode="lines", name=labels[4]),
scatter(x=V, y=i[:,5], mode="lines", name=labels[5])],
Layout(xaxis_range=[0,30], yaxis_range=[0,5], title="PV module - I(V)", xaxis_title="V", yaxis_title="I"))
V_ = V.*ones(5)';
P = i .* V_;
plot([
scatter(x=V, y=P[:,1], mode="lines", name=labels[1]),
scatter(x=V, y=P[:,2], mode="lines", name=labels[2]),
scatter(x=V, y=P[:,3], mode="lines", name=labels[3]),
scatter(x=V, y=P[:,4], mode="lines", name=labels[4]),
scatter(x=V, y=P[:,5], mode="lines", name=labels[5])],
Layout(xaxis_range=[0,30], yaxis_range=[0,100], title="PV module - P(V)", xaxis_title="V", yaxis_title="P"))

Side note regarding the V(I) characteristics:

For the V(I) characteristic, the axes are swapped during plotting and the explicit function is not used. When plotting with the explicit function, there is a problem with the logarithm, which must not be less than zero. Therefore, the range of values must be very fine tuned. But since we can see in the upper two plots how steep the function is, this is very difficult.

plot([
scatter(x=i[:,1], y=V, mode="lines", name=labels[1]),
scatter(x=i[:,2], y=V, mode="lines", name=labels[2]),
scatter(x=i[:,3], y=V, mode="lines", name=labels[3]),
scatter(x=i[:,4], y=V, mode="lines", name=labels[4]),
scatter(x=i[:,5], y=V, mode="lines", name=labels[5])],
Layout(xaxis_range=[0,5], yaxis_range=[0,30], title="PV module - V(I)", xaxis_title="I", yaxis_title="V"))
plot([
scatter(x=i[:,1], y=P[:,1], mode="lines", name=labels[1]),
scatter(x=i[:,2], y=P[:,2], mode="lines", name=labels[2]),
scatter(x=i[:,3], y=P[:,3], mode="lines", name=labels[3]),
scatter(x=i[:,4], y=P[:,4], mode="lines", name=labels[4]),
scatter(x=i[:,5], y=P[:,5], mode="lines", name=labels[5])],
Layout(xaxis_range=[0,5], yaxis_range=[0,80], title="PV module - P(I)", xaxis_title="I", yaxis_title="P"))

Select Test point

Here we use the functions that are also part of ElectricGrid.jl to evaluate the characteristic curves at selected points.

v_test = 25
i_test = get_I(ArrayPV, v_test, 1000, 27)
2.436914259130086
i_test2 = 2.436914259130086
_ = get_V(ArrayPV, i_test2, 1000, 27)
v_test2 = ArrayPV.v_next
25.0
Plot([scatter(x=V, y=i[:,5], mode="lines", name=labels[5]),
      scatter(x=[v_test], y=[i_test], mode="markers", name="Test point",marker=attr(size=10))], 
Layout(xaxis_range=[0,30], yaxis_range=[0,5], title="PV module - I(V)", xaxis_title="V", yaxis_title="I"))
Plot([scatter(x=i[:,5], y=V, mode="lines", name=labels[5]),
      scatter(x=[i_test2], y=[v_test2], mode="markers", name="Test point",marker=attr(size=10))],
Layout(xaxis_range=[0,5], yaxis_range=[0,30], title="PV module - V(I)", xaxis_title="I", yaxis_title="V"))

Modelling batteries

The modeling of the battery follows that described in the implementation in this papper. The battery is a lithium-ion battery, which is shown in the following equivalent circuit.

Equivalent curcuit diagram of the battery

The most important components are the input resistance $R_0$ and the RC pairs, which reflect the individual cells within the battery. The number of the cell is indicated by the index $n$. all quantities of the battery are dependent on the state of charge (SOC) and the temperature of the battery. Since these are nonlinear characteristics, they must be recorded experimentally and stored in look-up tables. Here, the values from the Matlab implementation were also used (Source).

Now we create a battery. To do this, we first store the characteristic curves for the individual cell. If no specific values are entered here, the default values described above are used.

Cell = battery_module();
Battery = battery_block(battery_module=Cell, parallel=1, serial=1);

Now that the battery has been created, a charge and final charge cycle is to be plotted in the following. By default, the SOC is set to 0. The battery can be accessed via the update function. With this function, the charge and discharge current can be defined and the temperature of the cell can be adjusted. The update function is called in every step.

I_s = 1; # Define a current of 13 A for charging and discharging

SOC_load = []
Voltage_load = []
Current_load = []

SOC_discharge = []
Voltage_discharge = []
Current_discharge = []
Any[]

The following cell might take some time, depending on your machine. To obtain a full charge and discharge cycle and plot the SOC requires a longer simulation, which will only be hinted at here:

# for i in 1:1e7 # Ts*Steps=T_total // 1e-4 * 1e7 * 7,4 =  
#     I_in = -I_s

#     get_V(Battery, I_in, 20)

#     append!(SOC_load, Battery.SOC)
#     append!(Current_load, I_in)
#     append!(Voltage_load, Battery.v_dc)
# end

# for i in 1:1e7
#     I_in = I_s

#     get_V(Battery, I_in, 20)

#     append!(SOC_discharge, Battery.SOC)
#     append!(Current_discharge, I_in)
#     append!(Voltage_discharge, Battery.v_dc)
# end

# plot([scatter(x=SOC_load[1:5000:end], y=Voltage_load[1:5000:end], mode="markers", name="Charge"),
# scatter(x=SOC_discharge[1:5000:end], y=Voltage_discharge[1:5000:end], mode="markers", name="Discharge")],
# Layout(title="Charge and discharge characteristic", xaxis_title="SOC", yaxis_title="V"))

Access the models

As already shown in the basic enviroment notebooks, we now want to create an enviroment that uses the presented models instead of the ideal sources. The first step is to create an env with the corresponding parameter dict entries, specifying the DC link. In the first example we want to have a fully charged battery being discharged over a load.

CM = [ 0. 1.
      -1. 0.]

R_load, L_load, X, Z = ParallelLoadImpedance(1e3, .95, 230)


parameters = Dict{Any,Any}(
    "source" => Any[
        Dict{Any,Any}("source_type" => "battery", "pwr" => 15000.0, "fltr" => "L",
                      "L1" => 0.001, "C" => 2e-8, "L2" => 0.001, "R_C" => 0.01,
                      "R2" => 0.05, "Q" => 100 * 26 * 3600, "mode" => 1,
                      "parallel" => 48, "serial" => 160, "i_bat_limit" => 1e3,
                      "i_limit" => 1e3, "v_limit" => 1e3),
            ],
    "load" => Any[
        Dict{Any,Any}("impedance" => "RL", "R" => R_load, "L" => L_load, "i_limit" => 500,
                      "v_limit" => 1000)
    ],
    "cable" => Any[
        Dict{Any,Any}("Cb" => 4.0e-9, "Lb" => 0.000264, "Rb" => 0.002, "C" => 0.4e-2, "i_limit" => 1.0e13, "v_limit" => 1000, "len" => 1.0, "L" => 0.25e-5, "R" => 0.1e-1),
    ],
    "grid" => Dict{Any,Any}("fs" => 1e4, "phase" => 3, "v_rms" => 230, "f_grid" => 50, "ramp_end" => 0.9)
)


states = []
states_all = []
Charging = []
Discharging = []
V_dc_char = []
V_dc_dis = []
env = ElectricGridEnv(num_sources=1, num_loads=1, CM=CM, parameters=parameters, t_end=2, verbosity=2)

Multi_Agent = SetupAgents(env)

hook = DataHook(collect_vdc_ids = [1],
                collect_soc_ids = [1],
                collect_idc_ids = [1],)

hook = Simulate(Multi_Agent, env, hook=hook);
[ Info: Normalization is done based on the defined parameter limits.
[ Info: Time simulation run time: 2.0 [s] ~> 20001 steps
[ Info: 1 'classically' controlled source has been initialised.
[ Info: 1 source has been set up in Swing mode.
[ Info: All 'classically' controlled sources have been automatically set up with droop coeficients, and proportional and integral gains.

Below we plot the results:

p = RenderHookResults(hook = hook,
                    vdc_to_plot     = [1],
                    soc_to_plot     = [1],
                    idc_to_plot     = [1],
                    return_plot     = false)

The struct vdc_link_voltages collects all information of the sources in a list:

env.vdc_link_voltages.sources
1-element Vector{Any}:
 BatteryVoltage(1, [1, 6, 11], [2, 7, 12], battery_block(battery_module([0.010342 0.0012244; 0.0067316 0.0011396; … ; 0.0032886 0.0010986; 0.0028114 0.0010309], [2287.7 11897.0; 6122.0 24515.0; … ; 12165.0 26430.0; 9118.0 24795.0]), 273.15, 0.016, [3.5042 3.5136; 3.5573 3.5646; … ; 4.0764 4.0821; 4.1924 4.193], 93593.98693730272, 0.9999357578771658, [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0], [293.15, 313.15], 0.0001, "discharge", 11×2 extrapolate(interpolate((::Vector{Float64},::Vector{Float64}), ::Matrix{Float64}, Gridded(Linear())), Throw()) with element type Float64:
 3.5042  3.5136
 3.5573  3.5646
 3.6009  3.6153
 3.6393  3.6565
 3.6742  3.6889
 3.7121  3.7214
 3.7937  3.8078
 3.8753  3.8945
 3.97    3.9859
 4.0764  4.0821
 4.1924  4.193, 11×2 extrapolate(interpolate((::Vector{Float64},::Vector{Float64}), ::Matrix{Float64}, Gridded(Linear())), Throw()) with element type Float64:
 0.010342   0.0012244
 0.0067316  0.0011396
 0.0051156  0.0012661
 0.0043447  0.0012265
 0.0038826  0.0011163
 0.0034226  0.0009968
 0.003346   0.0011458
 0.0033222  0.001345
 0.0033201  0.0013091
 0.0032886  0.0010986
 0.0028114  0.0010309, 1000.0, 17439.330026094067, 48, 160, 93600, 659.6212913723982, 659.6213009379546, 3.7629751757887195))

In the second example we want to create another PV array, which should supply a load.

CM = [ 0.  1.
      -1.  0.]

R_load, L_load, X, Z = ParallelLoadImpedance(1e3, 1., 230)

parameters = Dict{Any, Any}(
        "source" => Any[
                        Dict{Any, Any}("source_type" => "pv",
                        "fltr"=>"L", "L1"=>0.001, "C"=>2e-8, "L2"=>0.001, "R_C"=>0.01,
                        "R2"=> 0.05, "mode" => 1, "parallel" => 20, "serial" => 20,
                        "i_limit" => 1e4, "v_limit" => 1e5, "module_N_cell" => 36),
                        ],
        "load"   => Any[
                        Dict{Any, Any}("impedance" => "R", "R" => 1e100, "i_limit" => 10e4, "v_limit" => 10e4),
                        ],
        "cable"   => Any[
                        Dict{Any, Any}("R" => 1e-3, "L" => 1e-4, "C" => 1e-4, "i_limit" => 1e4, "v_limit" => 1e6,),
                        ],
        "grid" => Dict{Any, Any}("fs"=>1e4, "phase"=>3, "v_rms"=>230, "f_grid" => 50, "ramp_end" => 0.9)
    )


env = ElectricGridEnv(num_sources=1, num_loads=1, CM = CM, parameters = parameters, t_end=2, verbosity=2)

Multi_Agent = SetupAgents(env)

hook = DataHook(collect_state_ids = ["source1_i_L1_a", "source1_v_C_cables_a"],
                collect_vdc_ids=[1],
                collect_soc_ids=[1],
                collect_idc_ids=[1],)

hook = Simulate(Multi_Agent, env, hook=hook);
[ Info: Normalization is done based on the defined parameter limits.
[ Info: Time simulation run time: 2.0 [s] ~> 20001 steps
[ Info: 1 'classically' controlled source has been initialised.
[ Info: 1 source has been set up in Swing mode.
[ Info: All 'classically' controlled sources have been automatically set up with droop coeficients, and proportional and integral gains.
p = RenderHookResults(hook=hook,
                      states_to_plot  = ["source1_i_L1_a", "source1_v_C_cables_a"],
                      vdc_to_plot=[1],
                      idc_to_plot=[1],
                      return_plot=false);