IEEE39 Bus Tutorial - Part IV: Advanced Modeling & Parameter Optimization

This is the fourth and final part of the IEEE 39-bus tutorial series:

Part I: Model Creation - Build the network structure with buses, lines, and components
Part II: Initialization - Perform power flow calculations and dynamic initialization
Part III: Dynamic Simulation - Run time-domain simulations and analyze system behavior
Part IV: Advanced Modeling & Parameter Optimization (this tutorial) - Create custom components and optimize system parameters

In this tutorial, we'll demonstrate advanced PowerDynamics.jl capabilities by:

Creating a custom droop-controlled inverter component
Integrating it into the IEEE 39-bus system
Optimizing its parameters to improve system performance

This tutorial showcases custom component creation and the integration with Julia's optimization ecosystem for parameter tuning.

This script can be downloaded as a normal Julia script here.

Note

This tutorial is designed as a pedagogical example. It does not necessarily represent a realistic power system model and analysis, but rather serves to demonstrate the available tools while remaining relatively simple and concise.

# Loading required packages and setup
using PowerDynamics
using PowerDynamics.Library
using ModelingToolkit
using ModelingToolkit: D_nounits as Dt, t_nounits as t
using NetworkDynamics
using NetworkDynamics: SII
using OrdinaryDiffEqRosenbrock
using OrdinaryDiffEqNonlinearSolve
using SciMLSensitivity
using Optimization
using OptimizationOptimisers
using CairoMakie
using LinearAlgebra
using Graphs
using SparseConnectivityTracer
using Sparspak

# Load the network models from previous parts
EXAMPLEDIR = joinpath(pkgdir(PowerDynamics), "docs", "examples")
include(joinpath(EXAMPLEDIR, "ieee39_part1.jl"))  # Creates the basic network model
include(joinpath(EXAMPLEDIR, "ieee39_part3.jl"))  # Provides initialized network and reference solution

[ Info: Short circuit activated on line 5→8 at t = 0.1s
[ Info: Line 5→8 disconnected at t = 0.2s

Integration of a Droop-Controlled Inverter

In this section, we'll modify our network by adding a droop-controlled inverter.

Mathematical Background

The droop-controlled inverter implements a decentralized control strategy commonly used in microgrids and renewable energy integration. It establishes the following relationships:

Power Measurement:

\[\begin{aligned} P_{meas} &= u_r \cdot i_r + u_i \cdot i_i\\ Q_{meas} &= u_i \cdot i_r - u_r \cdot i_i \end{aligned}\]

Power Filtering (Low-pass filtering for measurement noise reduction):

\[\begin{aligned} \tau \cdot \frac{dP_{filt}}{dt} &= P_{meas} - P_{filt} \\ \tau \cdot \frac{dQ_{filt}}{dt} &= Q_{meas} - Q_{filt} \end{aligned}\]

Droop Control:

\[\begin{aligned} \omega &= \omega_0 - K_p \cdot (P_{filt} - P_{set}) \\ V &= V_{set} - K_q \cdot (Q_{filt} - Q_{set}) \end{aligned}\]

Voltage Angle Dynamics:

\[\frac{d\delta}{dt} = \omega - \omega_0\]

Output Voltage:

\[\begin{aligned} u_r &= V \cdot \cos(\delta) \\ u_i &= V \cdot \sin(\delta) \end{aligned}\]

These equations implement:

Frequency-Active Power Coupling (f-P): Frequency decreases when active power exceeds setpoint
Voltage-Reactive Power Coupling (V-Q): Voltage decreases when reactive power exceeds setpoint

This mimics the natural behavior of synchronous generators and enables stable power sharing in islanded operation.

Definition of the Droop Inverter Component

Network components in PowerDynamics must follow the Injector Interface - they connect to the network through a single Terminal:

      ┌──────────────────────────┐
(t)   │                          │
 o←───┤ Droop Inverter Equations │
      │                          │
      └──────────────────────────┘

We can define the following MTKModel to represent the droop inverter:

@mtkmodel DroopInverter begin
    @components begin
        terminal = Terminal()
    end

    @parameters begin
        Pset, [description="Active power setpoint", guess=1]
        Qset, [description="Reactive power setpoint", guess=0]
        Vset, [description="Voltage magnitude setpoint", guess=1]
        ω₀=1, [description="Nominal frequency"]
        Kp=1, [description="Active power droop coefficient"]
        Kq=0.1, [description="Reactive power droop coefficient"]
        τ_p = 1, [description="Active Power filter time constant"]
        τ_q = 1, [description="Reactive Power filter time constant"]
    end

    @variables begin
        Pmeas(t), [description="Measured active power", guess=1]
        Qmeas(t), [description="Measured reactive power", guess=0]
        Pfilt(t), [description="Filtered active power", guess=1]
        Qfilt(t), [description="Filtered reactive power", guess=1]
        ω(t), [description="Frequency"]
        δ(t), [description="Voltage angle", guess=0]
        V(t), [description="Voltage magnitude"]
    end

    @equations begin
        # Power measurement from terminal quantities
        Pmeas ~ terminal.u_r*terminal.i_r + terminal.u_i*terminal.i_i
        Qmeas ~ terminal.u_i*terminal.i_r - terminal.u_r*terminal.i_i

        # First-order low-pass filtering
        τ_p * Dt(Pfilt) ~ Pmeas - Pfilt
        τ_q * Dt(Qfilt) ~ Qmeas - Qfilt

        # Droop control equations
        ω ~ ω₀ - Kp * (Pfilt - Pset)  # Frequency decreases with excess power
        V ~ Vset - Kq * (Qfilt - Qset)  # Voltage decreases with excess reactive power

        # Voltage angle dynamics
        Dt(δ) ~ ω - ω₀

        # Output voltage components
        terminal.u_r ~ V*cos(δ)
        terminal.u_i ~ V*sin(δ)
    end
end;

Creating a Bus with the Droop Inverter

Following the descriptions in Modeling Concepts, we build an MTKBus using the droop as the single injector and then compile the bus model, similar to how we define the templates in part I:

           ╔═════════════════════════╗
           ║ Droop (compiled)        ║
 Network   ║  ┌────────────────────┐ ║
interface  ║  │ MTKBus             │ ║
 current ────→│┌──────┐ ┌────────┐ │ ║
           ║  ││BusBar├o┤Inverter│ │ ║
 voltage ←────│└──────┘ └────────┘ │ ║
           ║  └────────────────────┘ ║
           ╚═════════════════════════╝

@named inverter = DroopInverter()
mtkbus = MTKBus(inverter)
droop_bus_template = Bus(mtkbus; name=:DroopInverter)

VertexModel :DroopInverter NoFeedForward()
 ├─ 2 inputs:  [busbar₊i_r, busbar₊i_i]
 ├─ 3 states:  [inverter₊δ≈0, inverter₊Qfilt≈1, inverter₊Pfilt≈1]
 ├─ 2 outputs: [busbar₊u_r=1, busbar₊u_i=0]
 └─ 8 params:  [inverter₊Pset≈1, inverter₊Qset≈0, inverter₊Vset≈1, inverter₊ω₀=1, inverter₊Kp=1, inverter₊Kq=0.1, inverter₊τ_p=1, inverter₊τ_q=1]

We see that the droop inverter has 3 free parameters (you can check free_p(droop_bus_template) or dump_initial_state(droop_bus_template)). Therefore, similar to what we did in Part II, we need to help the initialization by attaching an additional initialization formula to the bus:

set_initformula!(
    droop_bus_template,
    @initformula(:inverter₊Vset = sqrt(:busbar₊u_r^2 + :busbar₊u_i^2))
)

Network Modification with Droop Inverter

We'll replace bus 32 (originally a controlled generator bus) with our new droop inverter bus.

DROOP_BUS_IDX = 32

To do so, we first collect all the "old" vertex and edge models.

Tip

We copy the components, to create individual instances for the new network. Since all the metadata (like the default values, initialize values and so on) are stored in the component models, we would otherwise share metadata between the old and the new network, which can lead to unexpected results.

vertex_models = [copy(nw[VIndex(i)]) for i in 1:nv(nw)];
edge_models = [copy(nw[EIndex(i)]) for i in 1:ne(nw)];

Now we need to replace the original bus model at index DROOP_BUS_IDX with our new droop inverter bus. However, we don't want to lose the original power flow model associated with this bus, so we need to attach it to the droop bus model:

original_pfmodel = get_pfmodel(vertex_models[DROOP_BUS_IDX])

We can then use the Bus constructor to essentially copy the droopbustemplate and adjust some properties, like the powerflow model and the vertex index.

droop_bus = Bus(droop_bus_template; pf=original_pfmodel, vidx=DROOP_BUS_IDX)

VertexModel :DroopInverter NoFeedForward() @ Vertex 32
 ├─ 2 inputs:  [busbar₊i_r, busbar₊i_i]
 ├─ 3 states:  [inverter₊δ≈0, inverter₊Qfilt≈1, inverter₊Pfilt≈1]
 ├─ 2 outputs: [busbar₊u_r=1, busbar₊u_i=0]
 ├─ 8 params:  [inverter₊Pset≈1, inverter₊Qset≈0, inverter₊Vset≈1, inverter₊ω₀=1, inverter₊Kp=1, inverter₊Kq=0.1, inverter₊τ_p=1, inverter₊τ_q=1]
 └─ 1 add. init eq. from 1 formula setting [:inverter₊Vset]
Powerflow model :pvbus with [pv₊P=6.5, pv₊V=0.9831]

We then replace the original bus model in the array with our droop bus and build a network again:

vertex_models[DROOP_BUS_IDX] = droop_bus
nw_droop = Network(vertex_models, edge_models)
set_jac_prototype!(nw_droop; remove_conditions=true)

Network with 189 states and 1272 parameters
 ├─ 39 vertices (6 unique types)
 └─ 46 edges (1 unique type)
Edge-Aggregation using SequentialAggregator(+)
Jacobian prototype defined: 3.45% sparsity
2 callback sets across 0 vertices and 1 edge

Additionally, we've set the jacobian prototype for performance gains during simulation and optimization, see NetworkDynamics docs on Sparsity Detection.

Network Initialization with Droop Inverter

The modified network requires the same initialization steps as the original:

Power flow solution
Dynamic component initialization

This all happens within initialize_from_pf!:

s0_droop = initialize_from_pf!(nw_droop; verbose=false)

Let's examine the initialized state of our droop inverter:

dump_initial_state(nw_droop[VIndex(DROOP_BUS_IDX)]; obs=false)

Inputs:
  busbar₊i_i            =  1.7883
  busbar₊i_r            = -6.6986
States:
  inverter₊Pfilt        =  6.5      (guess  1)
  inverter₊Qfilt        =  2.0514   (guess  1)
  inverter₊δ            =  0.044837 (guess  0)
Outputs:
  busbar₊u_i            =  0.044064
  busbar₊u_r            =  0.98211
Parameters:
  inverter₊Kp           =  1
  inverter₊Kq           =  0.1
  inverter₊Pset         =  6.5      (guess  1)
  inverter₊Qset         =  2.0514   (guess  0)
  inverter₊Vset         =  0.9831   (guess  1)
  inverter₊τ_p          =  1
  inverter₊τ_q          =  1
  inverter₊ω₀ (unused)  =  1

We see that the filtered powers match the setpoints (steady state), and both $P$ and $V_\mathrm{set}$ are initialized according to the parameters of the PV powerflow model.

Simulation with Droop Inverter

Now we'll simulate the modified network and compare it with the original system response:

prob_droop = ODEProblem(nw_droop, uflat(s0_droop), (0.0, 15.0), copy(pflat(s0_droop));
                       callback=get_callbacks(nw_droop))
sol_droop = solve(prob_droop, Rodas5P())
@assert SciMLBase.successful_retcode(sol_droop)

[ Info: Short circuit activated on line 5→8 at t = 0.1s
[ Info: Line 5→8 disconnected at t = 0.2s

Comparison of System Responses

Let's compare how the droop inverter affects the system's response to the short circuit disturbance:

let fig = Figure(; size=(1000, 600))
    selected_buses = [3, 4, 25, DROOP_BUS_IDX]  # Representative buses including the droop bus
    ts = range(0, 10, length=1000)

    for (i, bus) in enumerate(selected_buses)
        row, col = divrem(i-1, 2) .+ (1, 1)
        ax = Axis(fig[row, col];
                  title="Voltage Magnitude at Bus $bus" * (bus==DROOP_BUS_IDX ? " (droop bus)" : ""),
                  xlabel="Time [s]",
                  ylabel="Voltage [pu]")

        # Original system response
        lines!(ax, ts, sol(ts; idxs=VIndex(bus, :busbar₊u_mag)).u;
               label="Original System", color=:blue, linewidth=2)

        # Droop inverter system response
        lines!(ax, ts, sol_droop(ts; idxs=VIndex(bus, :busbar₊u_mag)).u;
               label="With Droop Inverter", color=:red, linewidth=2)

        ylims!(ax, 0.85, 1.15)
        i == 1 && axislegend(ax; position=:rb)
    end

    fig
end

We see that the overall system reacts similarly but distinctly differently to the identical disturbance.

Parameter Optimization

To showcase advanced capabilities of the SciML-ecosystem and the integration with PowerDynamics.jl, we now want to try to tune the droop inverter parameters so that the overall system behavior more closely resembles the original behavior, i.e., to reduce the difference between the system with generator and the system with droop inverter.

Optimization Problem Formulation

We define a loss function that measures the deviation between the original system response and the modified system response:

\[L(p) = \sum_{i,t} |x_{ref}(t)_i - x(t;p)_i|^2\]

Where we have

Parameters $p = [K_p, K_q, \tau]$ to be optimized
the reference solution $x_{ref}(t)$ (original system)
the solution of the modified system $x(t;p)$ with updated parameters $p$

Goal: Find parameters $p$ that minimize this loss function, making the droop inverter system behave as closely as possible to the original system.

Setting Up the Optimization

First, we define the reference solution and identify the tunable parameters.

We probe the original solution at fixed timepoints, exporting u_r and u_i for every bus:

opt_ref = sol(0.3:0.1:10, idxs=[VIndex(1:39, :busbar₊u_r), VIndex(1:39, :busbar₊u_i)])

Next, we need to identify the "tunable" parameters. This is a bit tricky, because the overall nw_droop has 1271 parameters, so we need to find the indices of the parameters we want to tune in the flat parameter array. We can do so, by leveraging NetworkDynamics implementation of the SymbolicIndexingInterface: tunable_parameters = [:inverter₊Kp, :inverter₊Kq, :inverter₊τ]

tunable_parameters = [:inverter₊Kp, :inverter₊Kq, :inverter₊τ_p, :inverter₊τ_q]
tp_idx = SII.parameter_index(sol_droop, VIndex(DROOP_BUS_IDX, tunable_parameters))

4-element Vector{Int64}:
 503
 504
 505
 506

We also get their initial values, which we use as the starting point for the optimization.

p0 = sol_droop(sol_droop.t[begin], idxs=collect(VIndex(DROOP_BUS_IDX, tunable_parameters)))

4-element Vector{Float64}:
 1.0
 0.1
 1.0
 1.0

Loss Function Implementation

The loss function simulates the system with given parameters and compares the result to the reference:

function loss(p)
    # Create parameter vector for the full system
    # allp = similar(p, pdim(nw_droop)) # create a vector of the full length
    allp = similar(p, length(s0_droop.p)) # create a vector of the full length
    allp .= pflat(s0_droop.p) # copy all "initial" parameters to that vector
    allp[tp_idx] .= p  # Update only the tunable parameters with the parameters for the given optimization iteration

    # Solve the system with new parameters
    _sol = solve(prob_droop, Rodas5P(autodiff=true);
        p = allp,
        saveat = opt_ref.t,
        tspan=(0.0, opt_ref.t[end]),
        initializealg = SciMLBase.NoInit(),
        abstol=0.01,
        reltol=0.01
    )

    # Return infinite loss if simulation failed
    if !SciMLBase.successful_retcode(_sol)
        @warn "Retcode $(_sol.retcode) indicates a failed simulation, returning Inf loss"
        return Inf
    end

    # Extract solution at reference time points
    x = _sol(opt_ref.t; idxs=[VIndex(1:39, :busbar₊u_r), VIndex(1:39, :busbar₊u_i)])

    # Compute L2 norm of the difference
    res = opt_ref.u - x.u
    l2loss = sum(abs2, reduce(vcat, res))
end

Optimization Execution

We use the Optimization.jl ecosystem with the Adam optimizer:

# Create optimization function with automatic differentiation
optf = Optimization.OptimizationFunction((x, p) -> loss(x), Optimization.AutoForwardDiff())

To better monitor the optimization progress, we want to store the optimized parameters at every iteration of the optimizer. We can do so by defining a callback function for the optimizer:

optimization_states = Any[] # global variable to store the optimization parameters at each step
callback = function (state, l)
    push!(optimization_states, state)
    println("Iteration $(state.iter): loss = $l\t p = $(state.u)")
    return false  # Continue optimization
end

That callback will snapshot the current parameter values at every step of the gradient descent.

With that, we can run the optimization:

optprob = Optimization.OptimizationProblem(optf, p0; callback)

@time optsol = Optimization.solve(optprob, Optimisers.Adam(0.06), maxiters = 50)

println("\nOptimization completed!")
println("Initial parameters: ", p0)
println("Optimized parameters: ", optsol.u)
println("Initial loss: ", loss(p0))
println("Final loss: ", loss(optsol.u))

Iteration 1: loss = 452.08075188374306	 p = [1.0, 0.1, 1.0, 1.0]
Iteration 2: loss = 309.92291410422524	 p = [1.059999999999443, 0.15999999999971232, 0.9400000000065613, 0.9400000000258409]
Iteration 3: loss = 226.82839704753798	 p = [1.119429899936713, 0.21256940171440714, 0.879929525043861, 0.947868891690498]
Iteration 4: loss = 175.6490386400381	 p = [1.1780147841755706, 0.25777646240837937, 0.8197515909382997, 0.9831253769802571]
Iteration 5: loss = 128.0715840762257	 p = [1.2355547875002968, 0.2963606610820473, 0.7593807157961712, 1.0282550822294012]
Iteration 6: loss = 89.87416214839155	 p = [1.2916704173826798, 0.3294265263474166, 0.6987841384665605, 1.0782158940049515]
Iteration 7: loss = 59.39735846960658	 p = [1.3459519376840596, 0.3579616972355813, 0.6380836254467024, 1.1306462126082624]
Iteration 8: loss = 35.55759311311985	 p = [1.3979563845396574, 0.3827887307862372, 0.577566686825637, 1.1838819343881544]
Iteration 9: loss = 19.289563386842634	 p = [1.4472095904450395, 0.40457492739364087, 0.5176550372901897, 1.2364203744906717]
Iteration 10: loss = 9.48257151282349	 p = [1.4932296277996917, 0.4238524817089607, 0.45910385915662016, 1.286805866604395]
Iteration 11: loss = 6.306839966619383	 p = [1.5354844284343463, 0.4410253169524846, 0.4034533495179298, 1.3336213443379286]
Iteration 12: loss = 9.408687455725875	 p = [1.5733917326665372, 0.4563830976254733, 0.3538267708543347, 1.3755900975997502]
Iteration 13: loss = 16.757691071284928	 p = [1.6063789993541222, 0.47010235131565203, 0.31520091344667717, 1.4118237116343528]
Iteration 14: loss = 24.6684872694376	 p = [1.6340650960419851, 0.4822826496228751, 0.2918118709392958, 1.4420159040055198]
Iteration 15: loss = 30.72718859941585	 p = [1.6563620560765537, 0.4930092257761865, 0.2844544196683683, 1.4662888416760076]
Iteration 16: loss = 33.29000163531275	 p = [1.673434920254641, 0.5023901890777868, 0.2910750894694208, 1.4849613387007996]
Iteration 17: loss = 32.58296667712008	 p = [1.6856810513711677, 0.5105773032215005, 0.3086430800475795, 1.4984056630989742]
Iteration 18: loss = 29.820749861222602	 p = [1.6936402483221973, 0.5177506941019066, 0.33440694343962746, 1.5070068208376641]
Iteration 19: loss = 25.25340868865076	 p = [1.6979007291495571, 0.5240893570334189, 0.3661640809291872, 1.5111863030520813]
Iteration 20: loss = 20.472909808720505	 p = [1.6991165207788734, 0.5297598948434372, 0.40200025356975333, 1.5114530449969512]
Iteration 21: loss = 15.73193524524867	 p = [1.697920747409479, 0.5349021984712309, 0.4403412970982441, 1.5083639447327355]
Iteration 22: loss = 12.2966586730192	 p = [1.694963315344509, 0.5396144185355319, 0.4795965112736686, 1.5026307183264247]
Iteration 23: loss = 9.812179548601593	 p = [1.690793354311234, 0.5439597843340049, 0.5184806739690873, 1.494951447082502]
Iteration 24: loss = 8.286578259510932	 p = [1.6858990712726074, 0.5479706223597108, 0.5559196154654787, 1.4860553494558078]
Iteration 25: loss = 7.513453844845163	 p = [1.6806980497812893, 0.5516579010158731, 0.5910911709875286, 1.4766589122650224]
Iteration 26: loss = 7.4774976583320445	 p = [1.6755315270309594, 0.5550205511874995, 0.6234503606890127, 1.4674252697366363]
Iteration 27: loss = 7.831846828343267	 p = [1.6706693317312973, 0.5580509186006186, 0.6526304490554632, 1.4589417656787274]
Iteration 28: loss = 8.384364187287433	 p = [1.6663191104858621, 0.5607400605318784, 0.6784238079982192, 1.451703471276441]
Iteration 29: loss = 9.052317722118543	 p = [1.6626503850282859, 0.5630792752736747, 0.7007259886077732, 1.4461371457705752]
Iteration 30: loss = 9.731201198534881	 p = [1.6597874189089408, 0.5650630601000995, 0.7195287442916565, 1.4425745828331964]
Iteration 31: loss = 10.44871216996037	 p = [1.6578099101174284, 0.5666960008268557, 0.7349004640713098, 1.4412123656848994]
Iteration 32: loss = 10.928094787267742	 p = [1.6567671605805185, 0.567982432715229, 0.746976175859379, 1.4421949000077958]
Iteration 33: loss = 11.234273615563305	 p = [1.6566599354914118, 0.5689348849214989, 0.7559279937693106, 1.445527097016171]
Iteration 34: loss = 11.365834210763712	 p = [1.6574580686774634, 0.5695715646617047, 0.7619588297659639, 1.451123544575911]
Iteration 35: loss = 11.332204746206843	 p = [1.6591026965165399, 0.5699159165235375, 0.7652957493359626, 1.4588184597488363]
Iteration 36: loss = 11.15530814475672	 p = [1.6615103832029234, 0.5699957761617619, 0.7661844695598227, 1.4683807785601943]
Iteration 37: loss = 10.866571769302723	 p = [1.664577723577094, 0.5698352840241352, 0.7648836183278228, 1.4795668954912076]
Iteration 38: loss = 10.505743421535144	 p = [1.668185637358339, 0.5694688496258944, 0.7616619104964079, 1.4920567097431159]
Iteration 39: loss = 10.110508973915483	 p = [1.6722042331794564, 0.5689317324313989, 0.7567925724269225, 1.5055122082742325]
Iteration 40: loss = 9.711341153020445	 p = [1.6764974934440822, 0.5682587291043054, 0.7505490407687287, 1.5195917457141033]
Iteration 41: loss = 9.33492085125624	 p = [1.6809281657717012, 0.5674830013533813, 0.7432012172606474, 1.533963014841813]
Iteration 42: loss = 9.004307498657155	 p = [1.6853618237646482, 0.5666351475430521, 0.7350130334012197, 1.5483133494414856]
Iteration 43: loss = 8.687302893336975	 p = [1.6896699207620969, 0.565742514540913, 0.7262406264893616, 1.5623572349599137]
Iteration 44: loss = 8.511020008457718	 p = [1.6937095076851032, 0.5648319662995943, 0.7171493970332993, 1.575791744826297]
Iteration 45: loss = 8.40440321440504	 p = [1.6973701785113366, 0.5639227642109207, 0.7079754308908606, 1.5884019677985508]
Iteration 46: loss = 8.358933569640651	 p = [1.7005546761745443, 0.5630300760008589, 0.698937485668068, 1.600011980776492]
Iteration 47: loss = 8.361022980493697	 p = [1.70318088966795, 0.5621651350208464, 0.6902331270962487, 1.6104858768992925]
Iteration 48: loss = 8.393893758114817	 p = [1.70518329811229, 0.5613356024731473, 0.6820353263193368, 1.6197271406771983]
Iteration 49: loss = 8.439859947525964	 p = [1.7065140840599544, 0.5605460583777222, 0.674488865902272, 1.6276769264025368]
Iteration 50: loss = 8.483888853319366	 p = [1.7071438667826457, 0.5597985283307083, 0.6677067036649383, 1.6343116494147467]
Iteration 50: loss = 6.306839966619383	 p = [1.5354844284343463, 0.4410253169524846, 0.4034533495179298, 1.3336213443379286]
 68.975666 seconds (147.48 M allocations: 38.089 GiB, 6.85% gc time, 59.55% compilation time: <1% of which was recompilation)

Optimization completed!
Initial parameters: [1.0, 0.1, 1.0, 1.0]
Optimized parameters: [1.5354844284343463, 0.4410253169524846, 0.4034533495179298, 1.3336213443379286]
Initial loss: 455.72825980863786
Final loss: 6.273945643671974

Optimization Results Analysis

Let's visualize how the optimization improved the system behavior:

function plot_optimization_comparison(p_initial, p_current)
    fig = Figure(; size=(1200, 800))
    selected_buses = [3, 4, 25, DROOP_BUS_IDX]
    ts = range(0, 10, length=1000)

    # Simulate with optimized parameters
    allp_opt = @lift let
        _p = copy(pflat(s0_droop))
        _p[tp_idx] .= $p_current
        _p
    end
    sol_opt = @lift solve(prob_droop, Rodas5P(); p=$allp_opt)


    for (i, bus) in enumerate(selected_buses)
        row, col = divrem(i-1, 2) .+ (1, 1)
        ax = Axis(fig[row, col];
                  title="Voltage Magnitude at Bus $bus" * (bus==DROOP_BUS_IDX ? " (droop bus)" : ""),
                  xlabel="Time [s]",
                  ylabel="Voltage [pu]")

        # Reference (original system)
        lines!(ax, ts, sol(ts; idxs=VIndex(bus, :busbar₊u_mag)).u;
               label="Reference", linestyle=:solid, color=:blue, linewidth=2)

        # Initial droop parameters
        lines!(ax, ts, sol_droop(ts; idxs=VIndex(bus, :busbar₊u_mag)).u;
               label="Initial Droop", linestyle=:dash, color=:red, linewidth=2)

        # Optimized droop parameters
        dat = @lift $(sol_opt)(ts; idxs=VIndex(bus, :busbar₊u_mag)).u
        lines!(ax, ts, dat; label="Optimized Droop", color=:green, linewidth=2)

        ylims!(ax, 0.85, 1.15)
        i == 1 && axislegend(ax; position=:rb)
    end

    fig
end

pobs = Observable(p0)
comparison_fig = plot_optimization_comparison(p0, pobs)
record(comparison_fig, "parameter_evolution.mp4", optimization_states; framerate=10) do s
    pobs[] = s.u
end

Parameter Evolution During Optimization

Let's see how each parameter changed during the optimization process:

let fig = Figure(; size=(1000, 800))
    param_names = ["Kp", "Kq", "τ_p", "τ_q"]
    for (i, param_name) in enumerate(param_names)
        row = (i-1) ÷ 2 + 1
        col = (i-1) % 2 + 1
        ax = Axis(fig[row, col];
                  title="Parameter Evolution: $param_name",
                  xlabel="Iteration",
                  ylabel="Parameter Value")

        # Extract parameter values over iterations
        param_values = [state.u[i] for state in optimization_states[1:end-1]]
        iterations = [state.iter for state in optimization_states[1:end-1]]

        scatterlines!(ax, iterations, param_values; linewidth=3, markersize=6, color=:blue)

        # Mark initial and final values
        hlines!(ax, [p0[i]]; linestyle=:dash, color=:gray, alpha=0.7)
        text!(ax, 1, p0[i]; text="Initial: $(round(p0[i], digits=3))",
              fontsize=10, color=:gray)

        lossax = Axis(fig[row, col],
            yticklabelcolor=:black,
            yaxisposition = :right,
            ylabel="loss", yscale=log10,
            xgridvisible=false, ygridvisible=false
        )
        scatterlines!(
            lossax,
            iterations,
            [loss(s.u) for s in optimization_states[1:end-1]],
            color=:black, linewidth=1, markersize=3,
        )
    end
    fig
end

This page was generated using Literate.jl.