diff --git a/build/ribasim_cli/Manifest.toml b/build/ribasim_cli/Manifest.toml index f87cbcb25..5fd921d52 100644 --- a/build/ribasim_cli/Manifest.toml +++ b/build/ribasim_cli/Manifest.toml @@ -886,7 +886,7 @@ uuid = "ae029012-a4dd-5104-9daa-d747884805df" version = "1.3.0" [[deps.Ribasim]] -deps = ["Arrow", "BasicModelInterface", "ComponentArrays", "Configurations", "DBInterface", "DataInterpolations", "DataStructures", "Dates", "Dictionaries", "DiffEqCallbacks", "Graphs", "IterTools", "Legolas", "OrdinaryDiffEq", "SQLite", "SciMLBase", "SparseArrays", "StructArrays", "Tables", "TimerOutputs"] +deps = ["Arrow", "BasicModelInterface", "ComponentArrays", "Configurations", "DBInterface", "DataInterpolations", "DataStructures", "Dates", "Dictionaries", "DiffEqCallbacks", "FiniteDiff", "Graphs", "IterTools", "Legolas", "Logging", "OrdinaryDiffEq", "SQLite", "SciMLBase", "SparseArrays", "StructArrays", "Tables", "TerminalLoggers", "TimerOutputs"] path = "../../core" uuid = "aac5e3d9-0b8f-4d4f-8241-b1a7a9632635" version = "0.1.0" diff --git a/core/Project.toml b/core/Project.toml index a4d7db22e..516b73bd7 100644 --- a/core/Project.toml +++ b/core/Project.toml @@ -14,15 +14,18 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e" Legolas = "741b9549-f6ed-4911-9fbf-4a1c0c97f0cd" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" SQLite = "0aa819cd-b072-5ff4-a722-6bc24af294d9" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" +TerminalLoggers = "5d786b92-1e48-4d6f-9151-6b4477ca9bed" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [compat] @@ -35,6 +38,7 @@ DataInterpolations = "3.7, 4" DataStructures = "0.18" Dictionaries = "0.3.25" DiffEqCallbacks = "2.23" +FiniteDiff = "2.21" Graphs = "1.6" IterTools = "1.4" Legolas = "0.5" @@ -43,6 +47,7 @@ SQLite = "1.5.1" SciMLBase = "1.60" StructArrays = "0.6.13" Tables = "1" +TerminalLoggers = "0.1.7" TimerOutputs = "0.5" julia = "1.9" diff --git a/core/src/Ribasim.jl b/core/src/Ribasim.jl index 851c7b330..d2deafd59 100644 --- a/core/src/Ribasim.jl +++ b/core/src/Ribasim.jl @@ -19,6 +19,7 @@ using SparseArrays using SQLite: SQLite, DB, Query, esc_id using StructArrays: StructVector using Tables: Tables, AbstractRow, columntable, getcolumn +using TerminalLoggers: TerminalLogger using TimerOutputs const to = TimerOutput() @@ -26,6 +27,7 @@ TimerOutputs.complement!() include("validation.jl") include("solve.jl") +include("jac.jl") include("config.jl") using .config: Config, Solver, algorithm, snake_case, zstd, lz4 include("utils.jl") diff --git a/core/src/bmi.jl b/core/src/bmi.jl index 6ac1ea4d6..dda0c0137 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -74,7 +74,7 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model timespan = (zero(t_end), t_end) jac_prototype = get_jac_prototype(parameters) - RHS = ODEFunction(water_balance!; jac_prototype) + RHS = ODEFunction(water_balance!; jac_prototype, jac = water_balance_jac!) @timeit_debug to "Setup ODEProblem" begin prob = ODEProblem(RHS, u0, timespan, parameters) @@ -209,13 +209,8 @@ function get_value( (; basin, flow_boundary) = p if variable == "level" - has_index, basin_idx = id_index(basin.node_id, feature_id) - - if !has_index - error("Level condition node #$feature_id is not a basin.") - end - - _, level = get_area_and_level(basin, basin_idx, storage[basin_idx]) + hasindex, basin_idx = id_index(basin.node_id, feature_id) + _, level, _ = get_area_and_level(basin, basin_idx, storage[basin_idx]) value = level elseif variable == "flow_rate" @@ -227,7 +222,7 @@ function get_value( value = flow_boundary.flow_rate[flow_boundary_idx](t) else - throw(ValueError("Unsupported condition variable $variable.")) + error("Unsupported condition variable $variable.") end return value diff --git a/core/src/create.jl b/core/src/create.jl index ec283e2f5..0d6b0be2d 100644 --- a/core/src/create.jl +++ b/core/src/create.jl @@ -298,6 +298,7 @@ function Basin(db::DB, config::Config)::Basin n = length(node_id) current_level = zeros(n) current_area = zeros(n) + current_darea = zeros(n) precipitation = fill(NaN, length(node_id)) potential_evaporation = fill(NaN, length(node_id)) @@ -326,6 +327,7 @@ function Basin(db::DB, config::Config)::Basin infiltration, current_level, current_area, + current_darea, area, level, storage, diff --git a/core/src/jac.jl b/core/src/jac.jl new file mode 100644 index 000000000..971836961 --- /dev/null +++ b/core/src/jac.jl @@ -0,0 +1,534 @@ +""" +The Jacobian is a n x n sparse matrix where n is the number of basins plus the number of +PidControl nodes. Each basin has a storage state, and each PidControl node has an error integral +state. If we write water_balance! as f(u, p(t), t) where u is the vector of all states, then +J[i,j] = ∂f_j/∂u_i. f_j dictates the time derivative of state j. + +J is very sparse because each state only depends on a small number of other states. +For more on the sparsity see get_jac_prototype. +""" +function water_balance_jac!( + J::SparseMatrixCSC{Float64, Int64}, + u::ComponentVector{Float64}, + p::Parameters, + t, +)::Nothing + (; basin) = p + J .= 0.0 + + # Ensures current_* vectors are current + set_current_basin_properties!(basin, u.storage, t) + + for nodefield in nodefields(p) + if nodefield == :pid_control + continue + end + + formulate_jac!(getfield(p, nodefield), J, u, p, t) + end + + # PID control must be done last + formulate_jac!(p.pid_control, J, u, p, t) + + return nothing +end + +""" +The contributions of LinearResistance nodes to the Jacobian. +""" +function formulate_jac!( + linear_resistance::LinearResistance, + J::SparseMatrixCSC{Float64, Int64}, + u::ComponentVector{Float64}, + p::Parameters, + t::Float64, +)::Nothing + (; basin, connectivity) = p + (; active, resistance, node_id) = linear_resistance + (; graph_flow) = connectivity + + for (id, isactive, R) in zip(node_id, active, resistance) + + # Inactive nodes do not contribute + if !isactive + continue + end + + id_in = only(inneighbors(graph_flow, id)) + id_out = only(outneighbors(graph_flow, id)) + + has_index_in, idx_in = id_index(basin.node_id, id_in) + has_index_out, idx_out = id_index(basin.node_id, id_out) + + if has_index_in + area_in = basin.current_area[idx_in] + term_in = 1 / (area_in * R) + J[idx_in, idx_in] -= term_in + end + + if has_index_out + area_out = basin.current_area[idx_out] + term_out = 1 / (area_out * R) + J[idx_out, idx_out] -= term_out + end + + if has_index_in && has_index_out + J[idx_in, idx_out] += term_out + J[idx_out, idx_in] += term_in + end + end + return nothing +end + +""" +The contributions of ManningResistance nodes to the Jacobian. +""" +function formulate_jac!( + manning_resistance::ManningResistance, + J::SparseMatrixCSC{Float64, Int64}, + u::ComponentVector{Float64}, + p::Parameters, + t::Float64, +)::Nothing + (; basin, connectivity) = p + (; node_id, active, length, manning_n, profile_width, profile_slope) = + manning_resistance + (; graph_flow) = connectivity + + for (i, id) in enumerate(node_id) + + # Inactive nodes do not contribute + if !active[i] + continue + end + + #TODO: this was copied from formulate! for the manning_resistance, + # maybe put in separate function + basin_a_id = only(inneighbors(graph_flow, id)) + basin_b_id = only(outneighbors(graph_flow, id)) + + h_a = get_level(p, basin_a_id) + h_b = get_level(p, basin_b_id) + bottom_a, bottom_b = basin_bottoms(basin, basin_a_id, basin_b_id, id) + slope = profile_slope[i] + width = profile_width[i] + n = manning_n[i] + L = length[i] + + Δh = h_a - h_b + q_sign = sign(Δh) + + # Average d, A, R + d_a = h_a - bottom_a + d_b = h_b - bottom_b + d = 0.5 * (d_a + d_b) + + A_a = width * d + slope * d_a^2 + A_b = width * d + slope * d_b^2 + A = 0.5 * (A_a + A_b) + + slope_unit_length = sqrt(slope^2 + 1.0) + P_a = width + 2.0 * d_a * slope_unit_length + P_b = width + 2.0 * d_b * slope_unit_length + R_h_a = A_a / P_a + R_h_b = A_b / P_b + R_h = 0.5 * (R_h_a + R_h_b) + + k = 1000.0 + kΔh = k * Δh + atankΔh = atan(k * Δh) + ΔhatankΔh = Δh * atankΔh + R_hpow = R_h^(2 / 3) + root = sqrt(2 / π * ΔhatankΔh) + + id_in = only(inneighbors(graph_flow, id)) + id_out = only(outneighbors(graph_flow, id)) + + has_index_in, idx_in = id_index(basin.node_id, id_in) + has_index_out, idx_out = id_index(basin.node_id, id_out) + + if has_index_in + basin_in_area = basin.current_area[idx_in] + ∂A_a = (width + 2 * slope * d_a) / basin_in_area + ∂A = 0.5 * ∂A_a + ∂P_a = 2 * slope_unit_length / basin_in_area + ∂R_h_a = (P_a * ∂A_a - A_a * ∂P_a) / P_a^2 + ∂R_h_b = width / (2 * basin_in_area * P_b) + ∂R_h = 0.5 * (∂R_h_a + ∂R_h_b) + # This float exact comparison is deliberate since `sqrt_contribution` has a + # removable singularity, i.e. it doesn't exist at $\Delta h = 0$ because of + # division by zero but the limit Δh → 0 does exist and is equal to the given value. + if Δh == 0 + sqrt_contribution = 2 / (sqrt(2 * π) * basin_in_area) + else + sqrt_contribution = + (atankΔh + kΔh / (1 + kΔh^2)) / + (basin_in_area * sqrt(2 * π * ΔhatankΔh)) + end + # term_in = q * (∂A / A + ∂R_h / R_h + sqrt_contribution) + term_in = + q_sign * ( + ∂A * R_hpow * root + + A * R_hpow * ∂R_h / R_h * root + + A * R_hpow * sqrt_contribution + ) / (n * sqrt(L)) + J[idx_in, idx_in] -= term_in + end + + if has_index_out + basin_out_area = basin.current_area[idx_out] + ∂A_b = (width + 2 * slope * d_b) / basin_out_area + ∂A = 0.5 * ∂A_b + ∂P_b = 2 * slope_unit_length / basin_out_area + ∂R_h_b = (P_b * ∂A_b - A_b * ∂P_b) / P_b^2 + ∂R_h_b = width / (2 * basin_out_area * P_b) + ∂R_h = 0.5 * (∂R_h_b + ∂R_h_a) + # This float exact comparison is deliberate since `sqrt_contribution` has a + # removable singularity, i.e. it doesn't exist at $\Delta h = 0$ because of + # division by zero but the limit Δh → 0 does exist and is equal to the given value. + if Δh == 0 + sqrt_contribution = 2 / (sqrt(2 * π) * basin_out_area) + else + sqrt_contribution = + (atankΔh + kΔh / (1 + kΔh^2)) / + (basin_out_area * sqrt(2 * π * ΔhatankΔh)) + end + # term_out = q * (∂A / A + ∂R_h / R_h + sqrt_contribution) + term_out = + q_sign * ( + ∂A * R_hpow * root + + A * R_hpow * ∂R_h / R_h * root + + A * R_hpow * sqrt_contribution + ) / (n * sqrt(L)) + J[idx_out, idx_out] -= term_out + end + + if has_index_in && has_index_out + J[idx_in, idx_out] += term_out + J[idx_out, idx_in] += term_in + end + end + return nothing +end + +""" +The contributions of Pump nodes to the Jacobian. +""" +function formulate_jac!( + pump::Pump, + J::SparseMatrixCSC{Float64, Int64}, + u::ComponentVector{Float64}, + p::Parameters, + t::Float64, +)::Nothing + (; basin, fractional_flow, connectivity, pid_control) = p + (; active, node_id, flow_rate, is_pid_controlled) = pump + + (; graph_flow) = connectivity + + for (i, id) in enumerate(node_id) + + # Inactive nodes do not contribute + if !active[i] + continue + end + + if is_pid_controlled[i] + continue + end + + id_in = only(inneighbors(graph_flow, id)) + + # For inneighbors only directly connected basins give a contribution + has_index_in, idx_in = id_index(basin.node_id, id_in) + + # For outneighbors there can be directly connected basins + # or basins connected via a fractional flow + # (but not both at the same time!) + if has_index_in + s = u.storage[idx_in] + + if s < 10.0 + dq = flow_rate[i] / 10.0 + + J[idx_in, idx_in] -= dq + + has_index_out, idx_out = id_index(basin.node_id, id_in) + + idxs_fractional_flow, idxs_out = get_fractional_flow_connected_basins( + id, + basin, + fractional_flow, + graph_flow, + ) + + # Assumes either one outneighbor basin or one or more + # outneighbor fractional flows + if isempty(idxs_out) + id_out = only(outneighbors(graph_flow, id)) + has_index_out, idx_out = id_index(basin.node_id, id_out) + + if has_index_out + J[idx_in, idx_out] += dq + end + else + for (idx_fractional_flow, idx_out) in + zip(idxs_fractional_flow, idxs_out) + J[idx_in, idx_out] += + dq * fractional_flow.fraction[idx_fractional_flow] + end + end + end + end + end + return nothing +end + +""" +The contributions of TabulatedRatingCurve nodes to the Jacobian. +""" +function formulate_jac!( + tabulated_rating_curve::TabulatedRatingCurve, + J::SparseMatrixCSC{Float64, Int64}, + u::ComponentVector{Float64}, + p::Parameters, + t::Float64, +)::Nothing + (; basin, fractional_flow, connectivity) = p + (; node_id, active, tables) = tabulated_rating_curve + (; graph_flow) = connectivity + + for (i, id) in enumerate(node_id) + if !active[i] + continue + end + + id_in = only(inneighbors(graph_flow, id)) + + # For inneighbors only directly connected basins give a contribution + has_index_in, idx_in = id_index(basin.node_id, id_in) + + # For outneighbors there can be directly connected basins + # or basins connected via a fractional flow + if has_index_in + # Computing this slope here is silly, + # should eventually be computed pre-simulation and cached! + table = tables[i] + levels = table.t + flows = table.u + level = basin.current_level[idx_in] + level_smaller_idx = searchsortedlast(table.t, level) + if level_smaller_idx == 0 + slope = 0.0 + else + if level_smaller_idx == length(flows) + level_smaller_idx = length(flows) - 1 + end + + slope = + (flows[level_smaller_idx + 1] - flows[level_smaller_idx]) / + (levels[level_smaller_idx + 1] - levels[level_smaller_idx]) + end + + dq = slope / basin.current_area[idx_in] + + J[idx_in, idx_in] -= dq + + idxs_fractional_flow, idxs_out = + get_fractional_flow_connected_basins(id, basin, fractional_flow, graph_flow) + + # Assumes either one outneighbor basin or one or more + # outneighbor fractional flows + if isempty(idxs_out) + id_out = only(outneighbors(graph_flow, id)) + has_index_out, idx_out = id_index(basin.node_id, id_out) + + if has_index_out + J[idx_in, idx_out] += dq + end + else + for (idx_fractional_flow, idx_out) in zip(idxs_fractional_flow, idxs_out) + J[idx_in, idx_out] += dq * fractional_flow.fraction[idx_fractional_flow] + end + end + end + end + return nothing +end + +""" +The contributions of TabulatedRatingCurve nodes to the Jacobian. +""" +function formulate_jac!( + pid_control::PidControl, + J::SparseMatrixCSC{Float64, Int64}, + u::ComponentVector{Float64}, + p::Parameters, + t::Float64, +)::Nothing + # TODO: update after pid_control equation test merge + (; basin, connectivity, pump) = p + (; node_id, active, listen_node_id, proportional, integral, derivative, error) = + pid_control + (; min_flow_rate, max_flow_rate) = pump + (; graph_flow, graph_control) = connectivity + + get_error!(pid_control, p) + + n_basins = length(basin.node_id) + integral_value = u.integral + + if any(.!isnan.(derivative)) + # Calling water_balance is expensive, but it is a sure way of getting + # the proper du for the pid controlled basins + # TODO: Do not allocate new memory here, make field of struct + du = zero(u) + water_balance!(du, u, p, t) + end + + for (i, id) in enumerate(node_id) + if !active[i] + continue + end + + id_pump = only(outneighbors(graph_control, id)) + listened_node_id = listen_node_id[i] + _, listened_node_idx = id_index(basin.node_id, listened_node_id) + listen_area = basin.current_area[listened_node_idx] + storage_listened_basin = u.storage[listened_node_idx] + area = basin.current_area[listened_node_idx] + reduction_factor = min(storage_listened_basin, 10.0) / 10.0 + + K_d = derivative[i] + if !isnan(K_d) + # dlevel/dstorage = 1/area + D = 1.0 - K_d * reduction_factor / area + else + D = 1.0 + end + + E = 0.0 + + K_p = proportional[i] + if !isnan(K_p) + E += K_p * error[i] + end + + K_i = integral[i] + if !isnan(K_i) + E += K_i * integral_value[i] + end + + if !isnan(K_d) + dtarget_level = 0.0 + du_listened_basin_old = du.storage[listened_node_idx] + E += K_d * (dtarget_level - du_listened_basin_old / area) + end + + # Clip values outside pump flow rate bounds + flow_rate = reduction_factor * E / D + was_clipped = false + + if flow_rate < min_flow_rate[i] + was_clipped = true + flow_rate = min_flow_rate[i] + end + + if !isnan(max_flow_rate[i]) + if flow_rate > max_flow_rate[i] + was_clipped = true + flow_rate = max_flow_rate[i] + end + end + + # PID control integral state + pid_state_idx = n_basins + i + J[pid_state_idx, listened_node_idx] -= 1 / listen_area + + # If the flow rate is clipped to one of the bounds it does + # not change with storages and thus doesn't contribute to the + # Jacobian + if was_clipped + continue + end + + # Only in this case the reduction factor has a non-zero derivative + reduction_factor_regime = (storage_listened_basin < 10.0) + + # Computing D and E derivatives + if !isnan(K_d) + darea = basin.current_darea[listened_node_idx] + dD = reduction_factor * darea + + if reduction_factor_regime + dD -= dreduction_factor / area + end + + dD *= K_d + + dE = + -K_d * ( + area * J[listened_node_idx, listened_node_idx] - + du.storage[listened_node_idx] * darea + ) / (area^2) + else + dD = 0.0 + dE = 0.0 + end + + if !isnan(K_p) + dE -= K_p / area + end + + if reduction_factor_regime + dreduction_factor = 0.1 + + dq = dreduction_factor * E / D + else + dq = 0.0 + end + + dq += reduction_factor * (D * dE - E * dD) / (D^2) + + J[listened_node_idx, listened_node_idx] -= dq + + id_out = only(outneighbors(graph_flow, id_pump)) + has_index, idx_out_out = id_index(basin.node_id, id_out) + + if has_index + J[pid_state_idx, idx_out_out] += K_i * reduction_factor / D + J[listened_node_idx, idx_out_out] += dq + end + end + return nothing +end + +""" +Method for nodes that do not contribute to the Jacobian +""" +function formulate_jac!( + node::AbstractParameterNode, + J::SparseMatrixCSC{Float64, Int64}, + u::ComponentVector{Float64}, + p::Parameters, + t::Float64, +)::Nothing + node_type = nameof(typeof(node)) + + if !isa( + node, + Union{ + Basin, + DiscreteControl, + FlowBoundary, + FractionalFlow, + LevelBoundary, + Terminal, + }, + ) + error( + "It is not specified how nodes of type $node_type contribute to the Jacobian.", + ) + end + return nothing +end diff --git a/core/src/solve.jl b/core/src/solve.jl index 245c6368c..4ebddea13 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -78,6 +78,9 @@ struct Basin{C} <: AbstractParameterNode # cache this to avoid recomputation current_level::Vector{Float64} current_area::Vector{Float64} + # The derivative of the area with respect to the level + # used for the analytical Jacobian + current_darea::Vector{Float64} # Discrete values for interpolation area::Vector{Vector{Float64}} level::Vector{Vector{Float64}} @@ -95,6 +98,7 @@ struct Basin{C} <: AbstractParameterNode infiltration, current_level, current_area, + current_darea, area, level, storage, @@ -111,6 +115,7 @@ struct Basin{C} <: AbstractParameterNode infiltration, current_level, current_area, + current_darea, area, level, storage, @@ -396,6 +401,20 @@ function valid_n_neighbors(graph::DiGraph{Int}, node::AbstractParameterNode)::Ve return errors end +function set_current_basin_properties!( + basin::Basin, + storage::AbstractVector{Float64}, + t::Real, +)::Nothing + for i in eachindex(storage) + s = storage[i] + area, level, darea = get_area_and_level(basin, i, s) + basin.current_level[i] = level + basin.current_area[i] = area + basin.current_darea[i] = darea + end +end + """ Linearize the evaporation flux when at small water depths Currently at less than 0.1 m. @@ -407,12 +426,11 @@ function formulate!( t::Real, )::Nothing for i in eachindex(storage) - s = storage[i] - area, level = get_area_and_level(basin, i, s) - basin.current_level[i] = level - basin.current_area[i] = area - bottom = basin.level[i][1] # add all precipitation that falls within the profile + level = basin.current_level[i] + area = basin.current_area[i] + + bottom = basin.level[i][1] fixed_area = basin.area[i][end] depth = max(level - bottom, 0.0) reduction_factor = min(depth, 0.1) / 0.1 @@ -482,36 +500,41 @@ function continuous_control!( flow_rate = 0.0 + K_d = derivative[i] + if !isnan(K_d) + # dlevel/dstorage = 1/area + area = basin.current_area[listened_node_idx] + D = 1.0 - K_d * reduction_factor / area + else + D = 1.0 + end + K_p = proportional[i] if !isnan(K_p) - flow_rate += reduction_factor * K_p * error[i] + flow_rate += reduction_factor * K_p * error[i] / D end K_i = integral[i] if !isnan(K_i) - flow_rate += reduction_factor * K_i * integral_value[i] + flow_rate += reduction_factor * K_i * integral_value[i] / D end - K_d = derivative[i] if !isnan(K_d) - # dlevel/dstorage = 1/area - area = basin.current_area[listened_node_idx] dtarget_level = 0.0 du_listened_basin_old = du.storage[listened_node_idx] # The expression below is the solution to an implicit equation for # du_listened_basin. This equation results from the fact that if the derivative # term in the PID controller is used, the controlled pump flow rate depends on itself. - du_listened_basin = - ( - du_listened_basin_old - flow_rate - - reduction_factor * K_d * dtarget_level - ) / (1 - reduction_factor * K_d / area) - flow_rate = du_listened_basin_old - du_listened_basin + flow_rate += K_d * (dtarget_level - du_listened_basin_old / area) / D end # Clip values outside pump flow rate bounds flow_rate = clamp(flow_rate, min_flow_rate[i], max_flow_rate[i]) + # Below du.storage is updated. This is normally only done + # in formulate!(du, connectivity, basin), but in this function + # flows are set so du has to be updated too. + pump.flow_rate[controlled_node_idx] = flow_rate du.storage[listened_node_idx] -= flow_rate @@ -522,10 +545,6 @@ function continuous_control!( flow[src_id, controlled_node_id] = flow_rate flow[controlled_node_id, dst_id] = flow_rate - # Below du.storage is updated. This is normally only done - # in formulate!(du, connectivity, basin), but in this function - # flows are set so du has to be updated too. - has_index, dst_idx = id_index(basin.node_id, dst_id) if has_index du.storage[dst_idx] += flow_rate @@ -632,7 +651,7 @@ The hydraulic radius is defined as: Where P is the wetted perimeter. -The "upstream" water depth is used to compute cross-sectional area and +The average of the upstream and downstream water depth is used to compute cross-sectional area and hydraulic radius. This ensures that a basin can receive water after it has gone dry. """ @@ -677,8 +696,9 @@ function formulate!(manning_resistance::ManningResistance, p::Parameters)::Nothi R_h_a = A_a / P_a R_h_b = A_b / P_b R_h = 0.5 * (R_h_a + R_h_b) + k = 1000.0 - q = q_sign * A / n * R_h^(2 / 3) * sqrt(Δh / L * 2 / π * atan(1000 * Δh)) + q = q_sign * A / n * R_h^(2 / 3) * sqrt(Δh / L * 2 / π * atan(k * Δh)) flow[basin_a_id, id] = q flow[id, basin_b_id] = q @@ -803,6 +823,9 @@ function formulate_flows!( return nothing end +""" +The right hand side function of the system of ODEs set up by Ribasim. +""" function water_balance!( du::ComponentVector{Float64}, u::ComponentVector{Float64}, @@ -817,7 +840,10 @@ function water_balance!( du .= 0.0 nonzeros(connectivity.flow) .= 0.0 - # ensures current_level is current + # Ensures current_* vectors are current + set_current_basin_properties!(basin, storage, t) + + # Basin forcings formulate!(du, basin, storage, t) # First formulate intermediate flows diff --git a/core/src/utils.jl b/core/src/utils.jl index 6e3686f59..2d502193c 100644 --- a/core/src/utils.jl +++ b/core/src/utils.jl @@ -55,11 +55,15 @@ function create_storage_tables( return area, level, storage end +""" +Compute the area and level of a basin given its storage. +Also returns darea/dlevel as it is needed for the Jacobian. +""" function get_area_and_level( basin::Basin, state_idx::Int, storage::Float64, -)::Tuple{Float64, Float64} +)::Tuple{Float64, Float64, Float64} storage_discrete = basin.storage[state_idx] area_discrete = basin.area[state_idx] level_discrete = basin.level[state_idx] @@ -72,7 +76,7 @@ function get_area_and_level( area_discrete::Vector{Float64}, level_discrete::Vector{Float64}, storage::Float64, -)::Tuple{Float64, Float64} +)::Tuple{Float64, Float64, Float64} # storage_idx: smallest index such that storage_discrete[storage_idx] >= storage storage_idx = searchsortedfirst(storage_discrete, storage) @@ -81,6 +85,13 @@ function get_area_and_level( level = level_discrete[1] area = area_discrete[1] + level_lower = level + level_higher = level_discrete[2] + area_lower = area + area_higher = area_discrete[2] + + darea = (area_higher - area_lower) / (level_higher - level_lower) + elseif storage_idx == length(storage_discrete) + 1 # With a storage above the profile, use a linear extrapolation of area(level) # based on the last 2 values. @@ -96,14 +107,14 @@ function get_area_and_level( if area_diff ≈ 0 # Constant area means linear interpolation of level + darea = 0.0 area = area_lower level = level_higher + level_diff * (storage - storage_higher) / (storage_higher - storage_lower) else - area = sqrt( - area_higher^2 + 2 * (storage - storage_higher) * area_diff / level_diff, - ) + darea = area_diff / level_diff + area = sqrt(area_higher^2 + 2 * (storage - storage_higher) * darea) level = level_lower + level_diff * (area - area_lower) / area_diff end @@ -120,19 +131,20 @@ function get_area_and_level( if area_diff ≈ 0 # Constant area means linear interpolation of level + darea = 0.0 area = area_lower level = level_lower + level_diff * (storage - storage_lower) / (storage_higher - storage_lower) else - area = - sqrt(area_lower^2 + 2 * (storage - storage_lower) * area_diff / level_diff) + darea = area_diff / level_diff + area = sqrt(area_lower^2 + 2 * (storage - storage_lower) * darea) level = level_lower + level_diff * (area - area_lower) / area_diff end end - return area, level + return area, level, darea end """ @@ -446,6 +458,7 @@ end """ If both the unique node upstream and the unique node downstream of these nodes are basins, then these directly depend on eachother and affect the Jacobian 2x +Basins always depend on themselves. """ function update_jac_prototype!( jac_prototype::SparseMatrixCSC{Float64, Int64}, @@ -456,15 +469,20 @@ function update_jac_prototype!( (; graph_flow) = connectivity for id in node.node_id - - # Only if the inneighbor and the outneighbor are basins - # do we get a contribution id_in = only(inneighbors(graph_flow, id)) id_out = only(outneighbors(graph_flow, id)) has_index_in, idx_in = id_index(basin.node_id, id_in) has_index_out, idx_out = id_index(basin.node_id, id_out) + if has_index_in + jac_prototype[idx_in, idx_in] = 1.0 + end + + if has_index_out + jac_prototype[idx_out, idx_out] = 1.0 + end + if has_index_in && has_index_out jac_prototype[idx_in, idx_out] = 1.0 jac_prototype[idx_out, idx_in] = 1.0 @@ -498,12 +516,14 @@ function update_jac_prototype!( "It is not specified how nodes of type $node_type contribute to the Jacobian prototype.", ) end + return nothing end """ If both the unique node upstream and the nodes down stream (or one node further if a fractional flow is in between) are basins, then the downstream basin depends on the upstream basin(s) and affect the Jacobian as many times as there are downstream basins +Upstream basins always depend on themselves. """ function update_jac_prototype!( jac_prototype::SparseMatrixCSC{Float64, Int64}, @@ -521,8 +541,11 @@ function update_jac_prototype!( # For outneighbors there can be directly connected basins # or basins connected via a fractional flow + # (but not both at the same time!) if has_index_in - idxs_out = + jac_prototype[idx_in, idx_in] = 1.0 + + _, idxs_out = get_fractional_flow_connected_basins(id, basin, fractional_flow, graph_flow) if isempty(idxs_out) @@ -532,10 +555,10 @@ function update_jac_prototype!( if has_index_out push!(idxs_out, idx_out) end - end - - for idx_out in idxs_out - jac_prototype[idx_in, idx_out] = 1.0 + else + for idx_out in idxs_out + jac_prototype[idx_in, idx_out] = 1.0 + end end end end @@ -587,14 +610,16 @@ function update_jac_prototype!( end """ -Get the state index of the basins that are connected to a node of given id via fractional flow. +Get the node type specific indices of the fractional flows and basins, +that are consecutively connected to a node of given id. """ function get_fractional_flow_connected_basins( node_id::Int, basin::Basin, fractional_flow::FractionalFlow, graph_flow::DiGraph{Int}, -)::Vector{Int} +)::Tuple{Vector{Int}, Vector{Int}} + fractional_flow_idxs = Int[] basin_idxs = Int[] for first_outneighbor_id in outneighbors(graph_flow, node_id) @@ -602,9 +627,13 @@ function get_fractional_flow_connected_basins( second_outneighbor_id = only(outneighbors(graph_flow, first_outneighbor_id)) has_index, basin_idx = id_index(basin.node_id, second_outneighbor_id) if has_index + push!( + fractional_flow_idxs, + searchsortedfirst(fractional_flow.node_id, first_outneighbor_id), + ) push!(basin_idxs, basin_idx) end end end - return basin_idxs + return fractional_flow_idxs, basin_idxs end diff --git a/core/test/bmi.jl b/core/test/bmi.jl index 03149ff30..1896c307f 100644 --- a/core/test/bmi.jl +++ b/core/test/bmi.jl @@ -26,7 +26,7 @@ end @testset "fixed timestepping" begin dict = to_dict(config_template) - dt = 3600 + dt = 10 dict["solver"] = Ribasim.Solver(; algorithm = "ImplicitEuler", adaptive = false, dt) config = from_dict(Ribasim.Config, dict) @test config.solver.algorithm == "ImplicitEuler" diff --git a/core/test/run_models.jl b/core/test/run_models.jl index d3f3a9072..a0296f053 100644 --- a/core/test/run_models.jl +++ b/core/test/run_models.jl @@ -48,14 +48,14 @@ end model = Ribasim.run(toml_path) @test model isa Ribasim.Model @test model.integrator.sol.retcode == Ribasim.ReturnCode.Success - @test model.integrator.sol.u[end] ≈ Float32[5.950373, 727.97125] skip = Sys.isapple() + @test model.integrator.sol.u[end] ≈ Float32[5.949285, 725.9446] skip = Sys.isapple() # the highest level in the dynamic table is updated to 1.2 from the callback @test model.integrator.p.tabulated_rating_curve.tables[end].t[end] == 1.2 end "Shorthand for Ribasim.get_area_and_level" function lookup(profile, S) - Ribasim.get_area_and_level(profile.S, profile.A, profile.h, S) + Ribasim.get_area_and_level(profile.S, profile.A, profile.h, S)[1:2] end @testset "Profile" begin @@ -172,5 +172,6 @@ end # https://www.hec.usace.army.mil/confluence/rasdocs/ras1dtechref/latest/theoretical-basis-for-one-dimensional-and-two-dimensional-hydrodynamic-calculations/1d-steady-flow-water-surface-profiles/friction-loss-evaluation @test all(isapprox.(h_expected, h_actual; atol = 0.02)) # Test for conservation of mass - @test all(isapprox.(model.saved_flow.saveval[end], 5.0)) skip = Sys.isapple() + @test all(isapprox.(model.saved_flow.saveval[end], 5.0, atol = 0.001)) skip = + Sys.isapple() end diff --git a/core/test/utils.jl b/core/test/utils.jl index 088680ec3..4492a5200 100644 --- a/core/test/utils.jl +++ b/core/test/utils.jl @@ -22,6 +22,7 @@ end # create two basins with different bottoms/levels area = [[0.01, 1.0], [0.01, 1.0]] level = [[0.0, 1.0], [4.0, 5.0]] + darea = zeros(2) storage = Ribasim.profile_storage.(level, area) target_level = [0.0, 0.0] basin = Ribasim.Basin( @@ -32,6 +33,7 @@ end [2.0, 3.0], [2.0, 3.0], [2.0, 3.0], + darea, area, level, storage, @@ -111,9 +113,9 @@ end @test jac_prototype.m == 4 @test jac_prototype.n == 4 - @test jac_prototype.colptr == [1, 2, 3, 4, 6] - @test jac_prototype.rowval == [2, 1, 2, 2, 3] - @test jac_prototype.nzval == ones(5) + @test jac_prototype.colptr == [1, 3, 5, 7, 9] + @test jac_prototype.rowval == [1, 2, 1, 2, 2, 3, 2, 4] + @test jac_prototype.nzval == ones(8) toml_path = normpath(@__DIR__, "../../data/pid_control/pid_control.toml") diff --git a/docs/contribute/addnode.qmd b/docs/contribute/addnode.qmd index 4dace81a5..59840749e 100644 --- a/docs/contribute/addnode.qmd +++ b/docs/contribute/addnode.qmd @@ -70,8 +70,21 @@ function formulate!(new_node_type::NewNodeType, p::Parameters)::Nothing end ``` -If the new node type can create dependendies between basins, this has to be taken into account in the sparsity structure of the Jacobian of `water_balance!` which is constructed in the function `get_jac_prototype` in `src\utils.jl`. +## The Jacobian +See [Equations](../core/equations.qmd#the-jacobian) for a mathematical description of the Jacobian. + +Before the Julia core runs its simulation, the sparsity structure `jac_prototype` of $J$ is determined with `get_jac_prototype` in `utils.jl`. This function runs trough all node types and looks for nodes that create dependencies between states. It creates a sparse matrix of zeros and ones, where the ones denote locations of possible non-zeros in $J$. + +We divide the various node types in groups based on what type of state dependencies they yield, and these groups are discussed below. Each group has its own method `update_jac_prototype!` in `utils.jl` for the sparsity structure induced by nodes of that group. `NewNodeType` should be added to the signature of one these methods, or to the list of node types that do not contribute to the Jacobian in the method of `update_jac_prototype!` whose signature contains `node::AbstractParameterNode`. Of course it is also possible that a new method of `update_jac_prototype!` has to be introduced. + +The current dependency groups are: + +- Out-neighbor dependencies: examples are `TabulatedRatingCurve`, `Pump` (the latter only in the reduction factor regime and not PID controlled). If the in-neighbor of a node of this group is a basin, then the storage of this basin affects itself and the storage of the outneighbor (or the basin one node further if it is connected with a `FractionalFlow` in between) if that is also a basin; +- Either-neighbor dependencies: examples are `LinearResistance`, `ManningResistance`. If either the in-neighbor or out-neighbor of a node of this group is a basin, the storage of this basin depends on itself. If both the in-neighbor and the out-neighbor are basins, their storages also depend on eachother. +- The `PidControl` node is a special case which is discussed in [equations](..\core\equations.qmd#sec-PID). + +In the methods `formulate_jac!` in `jac.jl` the analytical expressions for the partial derivatives $\frac{\partial Q_{i',j'}}{\partial u_i}$ (and the ones related to PID integral states) are hardcoded. For `NewNodeType` either a new method of `formulate_jac!` has to be introduced, or it has to be added to the list of node types that do not contribute to the Jacobian in the method of `formulate_jac!` whose signature contains `node::AbstractParameterNode`. # Python I/O diff --git a/docs/core/equations.qmd b/docs/core/equations.qmd index 4129fbfa2..8009dd789 100644 --- a/docs/core/equations.qmd +++ b/docs/core/equations.qmd @@ -43,14 +43,51 @@ Given a single basin with node ID $i \in B$, the equation that dictates the chan $$ \frac{\text{d}S_i}{\text{d}t} = -\sum_{(i',j') \in E | j' = i} Q_{i',j'} - \sum_{(i',j') \in E | i' = i} Q_{i',j'} + F_i(\mathbf{u},p,t). +\sum_{(i',j') \in E | j' = i} Q_{i',j'} - \sum_{(i',j') \in E | i' = i} Q_{i',j'} + F_i(p,t). $$ Here $Q_{i,j}$ is the flow along an edge, where the graph direction dictates positive flow. So the first term denotes flow towards the basin, the second one denotes flow away from the basin, and the third term denotes external forcing. -$F_i(\mathbf{u},p,t)$ is given by input data and the current state of the model, and $Q_{i' ,j'}$ is determined by the type of nodes that connect to that edge. +$F_i(p,t)$ is given by input data, and $Q_{i' ,j'}$ is determined by the type of nodes that connect to that edge. -The various node and forcing types that the model can contain are explained below. +The various node and forcing types that the model can contain are explained in the section [Natural water balance terms](equations.qmd#natural-water-balance-terms). + +::: {.callout-note} +In general a model has more nodes than states, so in the Julia core there is a distinction between node indices and state indices. For simplicity these are treated as equal in the documentation when it comes to basins and their storage. +::: + +## The Jacobian + +The Jacobian is a $n\times n$ matrix where $n$ is the number of states in the simulation. The computation of the Jacobian for the solver is hardcoded in the Julia core, which is faster and more accurate than using finite difference methods. + + +The entries of the Jacobian $J$ are given by +$$ +J[i,j] = \frac{\partial f_j}{\partial u_i}, +$$ + +hence $J[i,j]$ quantifies how $f_j$, the derivative of state $j$ with respect to time, changes with a change in state $i$. If a node creates dependendies between basin storages (or other states), then this yields contributions to the Jacobian. If $j$ corresponds to a storage state, then + +$$ +J[i,j] = \sum_{(i',j') \in E | j' = i} \frac{\partial Q_{i',j'}}{\partial u_i} - \sum_{(i',j') \in E | i' = i} \frac{\partial Q_{i',j'}}{\partial u_i}, +$$ + +Most of these terms are always $0$, because a flow over an edge only depends on a small number of states. Therefore the matrix $J$ is very sparse. + +For many contributions to the Jacobian the derivative of the level $l(S)$ of a basin with respect to its storage $S$ is required. To get an expression for this, we first look at the storage as a function of the level: + +$$ +S(l) = \int_{l_0}^l A(\ell)d\ell. +$$ + +From this we obtain $S'(l) = A(l)$ and thus +$$ +\frac{\text{d}l}{\text{d}S} = \frac{1}{A(S)}. +$$ + +:::{.callout-note} +The presence of division by the basin area means that areas of size zero are not allowed. +::: # Natural water balance terms ## Precipitation @@ -129,9 +166,11 @@ Ribasim's basins can be connected to each other, and each basin expects an explicit connection. These connections are currently available for inter-basin flows: -1. `OutflowTable` -2. `LinearResistance` -3. `ManningResistance` + +1. `Pump` +2. `TabulatedRatingCurve` +3. `LinearResistance` +4. `ManningResistance` The flow direction of the basin is not pre-determined: flow directions may freely reverse, provided the connection allows it. Currently, a `LinearResistance` @@ -144,19 +183,41 @@ basins (external nodes) in a network. 2. `LevelBoundary` 3. `FlowBoundary` -### OutflowTable +### Pump + +The behaviour of the pump is very straight forward if it is not PID controlled. Its flow is given by a fixed flow rate $q$, multiplied by a reduction factor: +$$ +Q_\text{pump} = \phi(u)q +$$ + +Here $u$ is the storage of the upstream basin. The reduction factor -The OutflowTable is a tabulation of a basin's discharge behavior. It describes -a piecewise linear relationship between the basin's storage volume and its +$$ +\phi(u) = +\begin{align} + \begin{cases} + \frac{u}{10} &\text{if}& 0 \leq u \le 10 \\ + 1 &\text{if}& u \geq 10 + \end{cases} +\end{align} +$$ {#eq-pumpreductionfactor} + +makes sure that the flow of the pump goes to $0$ as the upstream basin dries out. + +### TabulatedRatingCurve + +The Tabulated Rating Curve is a tabulation of a basin's discharge behavior. It describes +a piecewise linear relationship between the basin's level and its discharge. It can be understood as an empirical description of a basin's properties. This can include a weir, but also the lumped hydraulic behavior of the upstream channels. -The OutflowTable should indicate at which volume no discharge occurs (the dead + +The Tabulated Rating Curve should indicate at which volume no discharge occurs (the dead storage volume). :::{.callout-note} -Currently, the discharge relies only on the basin's volume; it could also use +Currently, the discharge relies only on the basin's level; it could also use the volume of both connected basins to simulate backwater effects, submersion of weirs, or even reversal of flows for high precipitation events. ::: @@ -342,38 +403,106 @@ $$ # PID controller {#sec-PID} -The PID controller continuously sets the flow rate of a pump or weir to bring the level of a certain basin closer to its setpoint. If we denote the setpoint by $r(t)$ and the basin level by $y(t)$, then the error is given by +The PID controller continuously sets the flow rate of a pump or weir to bring the level of a certain basin closer to its setpoint. If we denote the setpoint by $\text{SP}(t)$ and the basin level by $y(t)$, then the error is given by $$ -e(t) = r(t) - y(t). +e(t) = \text{SP}(t) - y(t). $$ {#eq-error} The output of the PID controller for the flow rate of the pump or weir is then given by $$ Q_\text{PID}(t) = K_p e(t) + K_i\int_{t_0}^t e(\tau)\text{d}\tau + K_d \frac{\text{d}e}{\text{d}t}, -$$ +$$ {#eq-PIDflow} for given constant parameters $K_p,K_i,K_d$. The pump or weir can have associated minimum and maximum flow rates $Q_\min, Q_\max$, and so $$ -Q_\text{pump/weir} = \begin{align} -\begin{cases} - Q_\min & \text{if} & Q_\text{PID} < Q_\min, \\ - Q_\text{PID} & \text{if} & Q_\min \leq Q_\text{PID} \leq Q_\max, \\ - Q_\max & \text{if} & Q_\text{PID} > Q_\max. -\end{cases} -\end{align} +Q_\text{pump/weir} = \text{clip}(\phi(u_\text{PID})Q_\text{PID}; Q_\min, Q_\max) $$ +where $u_\text{PID}$ is the storage of the controlled basin, $\phi$ is the reduction factor (see @eq-pumpreductionfactor) and + +$$ +\text{clip}(Q; Q_\min, Q_\max) = +\begin{align} + \begin{cases} + Q_\min & \text{if} & Q < Q_\min, \\ + Q & \text{if} & Q_\min \leq Q \leq Q_\max, \\ + Q_\max & \text{if} & Q > Q_\max. + \end{cases} +\end{align} +$$ For the integral term we denote $$ -I(t) = \int_{t_0}^t e(\tau)\text{d}\tau. +I(t) = \int_{t_0}^t e(\tau)\text{d}\tau, $$ -$I(t)$ is treated as a state of the system and thus it has its own equation in the system in @eq-system: +where $t_0$ is the last time the PID controller was made active. $I(t)$ is treated as a state of the system and thus it has its own equation in the system in @eq-system: $$ \frac{\text{d}I}{\text{d}t} = e(t). $$ +## The derivative term + +When $K_d \ne 0$ this adds a level of complexity. We can see this by looking at the error derivative more closely: +$$ +\frac{\text{d}e}{\text{d}t} = \frac{\text{d}\text{SP}}{\text{d}t} - \frac{1}{A(u_\text{PID})}\frac{\text{d}u_\text{PID}}{\text{d}t}, +$$ +where $A(u_\text{PID})$ is the area of the controlled basin. The complexity arises from the fact that $Q_\text{PID}$ is a contribution to $\frac{\text{d}u_\text{PID}}{\text{d}t} = f_\text{PID}$, which makes @eq-PIDflow an implicit equation for $Q_\text{PID}$. We define + +$$ +f_\text{PID} = \hat{f}_\text{PID} - Q_\text{pump/weir}, +$$ + +that is, $\hat{f}_\text{PID}$ is the right hand side of the ODE for the controlled basin storage state without the contribution of the PID controlled pump. The minus sign comes from the fact that the pump is puming away from the controlled basin. + +Using this, solving @eq-PIDflow for $Q_\text{PID}$ yields +$$ +Q_\text{pump/weir} = \text{clip}\left(\phi(u_\text{PID})\frac{K_pe + K_iI + K_d \left(\frac{\text{d}\text{SP}}{\text{d}t}-\frac{\hat{f}_\text{PID}}{A(u_\text{PID})}\right)}{1-\phi(u_\text{PID})\frac{K_d}{A(u_\text{PID})}};Q_\min,Q_\max\right), +$$ +where the clipping is again done last. Note that to compute this, $\hat{f}_\text{PID}$ has to be known first, meaning that the PID controlled pump flow rate has to be computed after all other contributions to the PID controlled basin's storage are known. + +## Jacobian contributions + +For convenience here we denote a time derivative of a variable as a dot over te symbol. + +For the integral state we simply get the contribution +$$ +\frac{\partial \dot{I}}{\partial u_\text{PID}} = - \frac{1}{A(u_\text{PID})}. +$$ + +Note that when the calculated pump flow lies outside $[Q_\min, Q_\max]$, the pump flow is locally constant and thus all partial derivatives of $Q_\text{pump/weir}$ with respect to states are $0$. Assuming that it lies in $[Q_\min, Q_\max]$, we derive the Jacobian contributions. Define the enumerator and denominator of the fraction in $Q_\text{pump/weir}$: + +$$ +\begin{align} + E(u_\text{PID}) &= K_pe + K_iI + K_d \left(\dot{\text{SP}}-\frac{\hat{f}_\text{PID}}{A(u_\text{PID})}\right) \\ + D(u_\text{PID}) &= 1-\phi(u_\text{PID})\frac{K_d}{A(u_\text{PID})} +\end{align} +$$ + +Then + +$$ +\frac{\partial Q_\text{pump/weir}}{\partial I} = K_i\frac{\phi(u_\text{PID})}{D(u_\text{PID})}. +$$ + +For the derivative with respect to $u_\text{PID}$ we need the derivative of the enumerator and denominator: +$$ +\begin{align} + E'(u_\text{PID}) &= -\frac{K_p}{A(u_\text{PID})} - K_d\frac{A(u_\text{PID})\frac{\partial\hat{f}_\text{PID}}{\partial u_\text{PID}} - \hat{f}_\text{PID}A'(u_\text{PID})}{A(u_\text{PID})^2}\\ + D'(u_\text{PID}) &= - K_d\left[\frac{\phi'(u_\text{PID})}{A(u_\text{PID})}- \phi(u_\text{PID})\frac{A'(u_\text{PID})}{A(u_\text{PID})^2}\right] +\end{align} +$$ +For computational efficiency we exploit that $\phi'(u_\text{PID})$ is mostly $0$. Note that + +$$ + \frac{\partial\hat{f}_\text{PID}}{\partial u_\text{PID}} = \hat{J}[\text{PID},\text{PID}], +$$ +\emph{i.e.} the Jacobian term of the dependence of $u_\text{PID}$ on itself without the contribution of the PID controlled pump/weir. Note thus that also the Jacobian contribution of the PID controlled pump can only be computed after the other contributions to the Jacobian of $u_\text{PID}$ to itself are known. + +The partial derivative with respect to $u_\text{PID}$ is then given by +$$ +\frac{\partial Q_\text{pump/weir}}{\partial u_\text{PID}} = \phi'(u_\text{PID})\frac{E(u_\text{PID})}{D(u_\text{PID})} + \phi(u_\text{PID})\frac{D(u_\text{PID})E'(u_\text{PID}) - E(u_\text{PID})D'(u_\text{PID})}{D(u_\text{PID})^2}. +$$ ## The sign of the parameters diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd index aa754790e..6f5bb6fca 100644 --- a/docs/core/usage.qmd +++ b/docs/core/usage.qmd @@ -455,7 +455,7 @@ control_state | String ### PidControl -The PidControl node controls the level in a basin by continuously controlling the flow rate of a connected pump or weir. See also [PID controller](https://en.wikipedia.org/wiki/PID_controller). When A PidControl node is made inactive, the node under its control retains the last flow rate value, and the error integral is reset to 0. +The PidControl node controls the level in a basin by continuously controlling the flow rate of a connected pump or weir. See also [PID controller](https://en.wikipedia.org/wiki/PID_controller). When A PidControl node is made inactive, the node under its control retains the last flow rate value, and the error integral is reset to 0. For computational efficiency, if one of the terms should not be used, set its coefficient to `NaN` in stead of $0$. column | type | unit | restriction -------------- | -------- | -------- | ----------- diff --git a/docs/python/examples.ipynb b/docs/python/examples.ipynb index 956c8e557..8b75b689f 100644 --- a/docs/python/examples.ipynb +++ b/docs/python/examples.ipynb @@ -753,7 +753,7 @@ " }\n", ")\n", "\n", - "state = pd.DataFrame(data={\"node_id\": [1, 3], \"storage\": [100.0, 0.0]})\n", + "state = pd.DataFrame(data={\"node_id\": [1, 3], \"storage\": [100.0, 1e-3]})\n", "\n", "basin = ribasim.Basin(profile=profile, static=static, state=state)" ] @@ -998,7 +998,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.11.4" } }, "nbformat": 4, diff --git a/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py b/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py index a2224ea5a..5cc8581b0 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py +++ b/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py @@ -73,7 +73,7 @@ def pump_discrete_control_model() -> ribasim.Model: } ) - state = pd.DataFrame(data={"node_id": [1, 3], "storage": [100.0, 0.0]}) + state = pd.DataFrame(data={"node_id": [1, 3], "storage": [100.0, 1e-3]}) basin = ribasim.Basin(profile=profile, static=static, state=state) diff --git a/python/ribasim_testmodels/ribasim_testmodels/pid_control.py b/python/ribasim_testmodels/ribasim_testmodels/pid_control.py index 1bed81234..ad118c8ec 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/pid_control.py +++ b/python/ribasim_testmodels/ribasim_testmodels/pid_control.py @@ -134,7 +134,7 @@ def pid_control_model(): pump=pump, pid_control=pid_control, starttime="2020-01-01 00:00:00", - endtime="2021-01-01 00:00:00", + endtime="2020-07-01 00:00:00", ) return model