From 6b969ac1f0ae7d75e20124a9ee325b69190c6fbd Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Tue, 13 Aug 2024 18:07:19 -0700 Subject: [PATCH] Adding SOS state vector initialization in gas inj While computing gas injection, when low pass filter is applied to the flow rate, it would show a transient if the flow rate starts at a non-zero value. This could be bad for MPC which uses the computation while optimization. To avoid this, we initialize the state vector of the low pass filter to the initial value of the flow rate. This way, there will never be an initial transient due to the low pass filter. --- src/gas_injection.jl | 32 +++++++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/src/gas_injection.jl b/src/gas_injection.jl index 4bc2de2..91619bd 100644 --- a/src/gas_injection.jl +++ b/src/gas_injection.jl @@ -230,7 +230,7 @@ function compute_gas_injection( ) end if !isnothing(LPF) - flow_rate = filt(LPF, flow_rate) + flow_rate = filt(LPF, flow_rate, create_si(LPF, flow_rate[1])) end flow_rate = map((x)::Float64 -> x < 0.0 ? 0.0 : x, flow_rate) future_flow_rates = flow_rate[end-skip+2:end] @@ -345,6 +345,36 @@ function get_lpf( return convert(SecondOrderSections, bilinear(filter_s, fs)) end +""" + create_si(filter, last_val) + +Create initial state for a filter based on a provided last value. This assumes that +upto this point, there was a constant command on (called prev_xi in this code) and the +filter had reached equilibrium to a constant output which is the last_val provided. +""" +function create_si(filter::SecondOrderSections, last_val::Float64) + n = length(filter.biquads) + g = filter.g + si = zeros(Float64, 2, n) + for fi ∈ n:-1:1 + b0 = filter.biquads[fi].b0 + b1 = filter.biquads[fi].b1 + b2 = filter.biquads[fi].b2 + a1 = filter.biquads[fi].a1 + a2 = filter.biquads[fi].a2 + if fi == n + prev_yi = last_val / g + else + prev_yi = last_val + end + prev_xi = prev_yi * (1 + a1 + a2) / (b0 + b1 + b2) + si[2, fi] = b2 * prev_xi - a2 * prev_yi + si[1, fi] = si[2, fi] + b1 * prev_xi - a1 * prev_yi + last_val = prev_xi + end + return si +end + """ dribble( data::Vector{Float64},