From c86acaedd45ec3330034a2486b610d949dbf1fc6 Mon Sep 17 00:00:00 2001 From: Philip Top Date: Fri, 28 Aug 2026 17:00:08 -0700 Subject: [PATCH 1/6] update the block library, fix some potential issues and clean up implementation and documentation --- src/griddyn/Block.h | 23 +- src/griddyn/CMakeLists.txt | 1 + src/griddyn/blocks/Block.cpp | 6 + src/griddyn/blocks/BlockSequence.h | 17 +- src/griddyn/blocks/ControlBlock.h | 17 +- src/griddyn/blocks/DeadbandBlock.h | 19 +- src/griddyn/blocks/DelayBlock.h | 56 +- src/griddyn/blocks/DerivativeBlock.h | 14 +- src/griddyn/blocks/FilteredDerivativeBlock.h | 14 +- src/griddyn/blocks/FunctionBlock.cpp | 5 +- src/griddyn/blocks/FunctionBlock.h | 14 +- src/griddyn/blocks/IntegralBlock.h | 13 +- src/griddyn/blocks/LeadLag.h | 124 ++++ src/griddyn/blocks/LeadLagBlock.cpp | 206 +++++++ src/griddyn/blocks/LeadLagBlock.h | 94 ++++ src/griddyn/blocks/LutBlock.cpp | 201 +++++-- src/griddyn/blocks/LutBlock.h | 42 +- src/griddyn/blocks/NullBlock.h | 8 +- src/griddyn/blocks/PidBlock.h | 14 +- src/griddyn/blocks/RampLimiter.h | 13 +- src/griddyn/blocks/TransferFunctionBlock.cpp | 561 +++++++++++++------ src/griddyn/blocks/TransferFunctionBlock.h | 60 +- src/griddyn/blocks/ValueLimiter.h | 13 +- src/griddyn/blocks/blockLibrary.h | 1 + src/griddyn/governors/GovernorTgov1.cpp | 35 +- src/griddyn/governors/GovernorTgov1.h | 15 +- test/componentTests/testBlocks.cpp | 128 ++++- 27 files changed, 1374 insertions(+), 340 deletions(-) create mode 100644 src/griddyn/blocks/LeadLag.h create mode 100644 src/griddyn/blocks/LeadLagBlock.cpp create mode 100644 src/griddyn/blocks/LeadLagBlock.h diff --git a/src/griddyn/Block.h b/src/griddyn/Block.h index cfedfc15d..c2767411d 100644 --- a/src/griddyn/Block.h +++ b/src/griddyn/Block.h @@ -17,15 +17,20 @@ namespace blocks { class RampLimiter; } // namespace blocks -/** @brief class implementing basic control system block - the basic block class takes a single input X the output is then \f$K*(X+bias)\f$ -optionally implementing limiters Omax and Omin the limiters have a reset level specified by -resetLevel once the object is initialized the determination of whether to use the ramps is fixed and -cannot be changed unless the object is reinitialized directly - -the blocks take 1 or 2 inputs the first being the single input, if the differential input is set -then the second argument is the time derivative of the input -*/ +/** + * @brief Common single-input/single-output base for dynamic control blocks. + * + * In its basic algebraic form GridBlock evaluates @f$y=K(u+b)@f$. Derived + * blocks select whether they own an algebraic or differential output/state, + * then override the residual, derivative, and analytic Jacobian as needed. + * The optional second input is @f$\dot u@f$ when `differential_input` is set. + * + * Inherited output limits apply a @ref blocks::ValueLimiter; inherited ramp + * limits apply a @ref blocks::RampLimiter to a differential output. Both use + * root-triggered engagement and `resetlevel` hysteresis. State layout and + * whether limits are active are fixed during dynamic initialization, so limit + * flags must be changed before reinitializing the block. + */ class GridBlock: public GridSubModel { public: /** @brief flags common for all control blocks diff --git a/src/griddyn/CMakeLists.txt b/src/griddyn/CMakeLists.txt index 3fd7acc6a..cbcb6abea 100644 --- a/src/griddyn/CMakeLists.txt +++ b/src/griddyn/CMakeLists.txt @@ -168,6 +168,7 @@ set(block_sources blocks/DerivativeBlock.cpp blocks/FunctionBlock.cpp blocks/LutBlock.cpp + blocks/LeadLagBlock.cpp blocks/TransferFunctionBlock.cpp blocks/FilteredDerivativeBlock.cpp blocks/BlockSequence.cpp diff --git a/src/griddyn/blocks/Block.cpp b/src/griddyn/blocks/Block.cpp index 1686e3c88..a795c26c4 100644 --- a/src/griddyn/blocks/Block.cpp +++ b/src/griddyn/blocks/Block.cpp @@ -22,6 +22,8 @@ static const TypeFactory BLOCK_FACTORY("block", std::to_array({"basic", "gain"}), "basic"); static const ChildTypeFactory CONTROL_BLOCK_FACTORY("block", "control"); +static const ChildTypeFactory + LEAD_LAG_BLOCK_FACTORY("block", std::to_array({"leadlag", "lead_lag"})); static const ChildTypeFactory DEADBAND_BLOCK_FACTORY("block", std::to_array({"deadband", "db"})); static const ChildTypeFactory @@ -40,6 +42,10 @@ static const ChildTypeFactory FILTERED_DERIVATIVE_BLOCK_FACTORY( "block", std::to_array({"fder", "filtered_deriv", "filtered_derivative"})); +static const ChildTypeFactory + TRANSFER_FUNCTION_BLOCK_FACTORY( + "block", + std::to_array({"transfer_function", "transferfunction", "tf"})); GridBlock::GridBlock(const std::string& objName): GridSubModel(objName) { diff --git a/src/griddyn/blocks/BlockSequence.h b/src/griddyn/blocks/BlockSequence.h index 5ed406c87..88f627e45 100644 --- a/src/griddyn/blocks/BlockSequence.h +++ b/src/griddyn/blocks/BlockSequence.h @@ -11,10 +11,19 @@ #include namespace griddyn::blocks { -/** @brief class implementing a sequence of blocks as a single block -A block is defined as a single input single output subModel. This object takes any number of blocks -in a sequence and processes them in the appropriate fashion. -*/ +/** + * @brief Ordered serial composition of one or more GridBlock instances. + * + * For blocks @f$B_0,\ldots,B_n@f$, this composite evaluates + * @f$y=B_n(\ldots B_1(B_0(u+b))\ldots)@f$. It propagates output derivatives, + * residuals, analytic Jacobians, roots, and root triggers through the same + * ordering. States remain owned by the child blocks; this class only supplies + * the single-input/single-output composition and local-step scheduling. + * + * `differential_input` causes the first child to receive an input derivative. + * BlockSequence is a composite container, not an independently factory-loaded + * elementary block. + */ class BlockSequence: public GridBlock { public: protected: diff --git a/src/griddyn/blocks/ControlBlock.h b/src/griddyn/blocks/ControlBlock.h index 3382869b4..79181e4eb 100644 --- a/src/griddyn/blocks/ControlBlock.h +++ b/src/griddyn/blocks/ControlBlock.h @@ -10,11 +10,18 @@ #include namespace griddyn::blocks { -/** @brief class implementing a control block -block implementing \f$H(S)=\frac{K(1+T_2 s}{1+T_1 s}\f$ -default is \f$T_2 =0\f$ for behavior equivalent to a delay block -if T1 is 0 it behaves like the basic block -*/ +/** + * @brief Legacy first-order lead--lag GridBlock. + * + * This predates @ref LeadLagBlock and implements + * @f$H(s)=K(1+T_2s)/(1+T_1s)@f$. It owns an intermediate differential + * state and exposes an algebraic output. `t1` is the denominator time + * constant and `t2` is the numerator time constant; `t2=0` gives the + * first-order lag form. Inherited GridBlock gain, bias, and limiter settings + * apply. New controller code should prefer @ref LeadLagBlock, whose equation, + * initialization contract, and zero-order-hold stepping are documented and + * tested explicitly. + */ class ControlBlock: public GridBlock { public: protected: diff --git a/src/griddyn/blocks/DeadbandBlock.h b/src/griddyn/blocks/DeadbandBlock.h index 57301f731..1077693df 100644 --- a/src/griddyn/blocks/DeadbandBlock.h +++ b/src/griddyn/blocks/DeadbandBlock.h @@ -11,9 +11,22 @@ #include namespace griddyn::blocks { -/** @brief class implementing a deadband system -TOBE added -*/ +/** + * @brief Stateful deadband, with optional continuous transition ramps. + * + * For biased input @f$x=u+b@f$, the normal deadband output is the configured + * `level` while @f$x\in[low,high]@f$ and is @f$x@f$ outside that interval. + * With the `shifted` flag, the outside branches are offset to meet `level` + * continuously at the bounds. `ramp`, `rampup`, and `rampdown` introduce + * linear transition regions instead. The exposed output is @f$K f(x)@f$ and + * inherited GridBlock limits may further clamp it. + * + * Root events transition among NORMAL, OUTSIDE, SHIFTED, and the two ramp + * states. `reset`, `resethigh`, and `resetlow` provide hysteresis; this is + * important for preventing event chatter near a boundary. `db`/`deadband` + * sets symmetric bounds around `level`; `high` and `low` set individual + * boundaries. + */ class DeadbandBlock: public GridBlock { public: /** @brief flags for the deadband block*/ diff --git a/src/griddyn/blocks/DelayBlock.h b/src/griddyn/blocks/DelayBlock.h index 7d34b225d..daaa5ec2c 100644 --- a/src/griddyn/blocks/DelayBlock.h +++ b/src/griddyn/blocks/DelayBlock.h @@ -6,31 +6,57 @@ #pragma once +/** + * @file DelayBlock.h + * @brief First-order lag (measurement or transport-delay approximation) block. + */ + #include "../Block.h" #include namespace griddyn::blocks { -/** @brief class implementing a delay block -block implementing \f$H(S)=\frac{K}{1+T_1 s}\f$ -if the time constant is very small it reverts to the basic block -*/ +/** + * @brief First-order lag block. + * + * This block realizes the proper transfer function + * @f[ + * H(s)=\frac{K}{1+T_1s}. + * @f] + * For input @f$u@f$, GridBlock bias @f$b@f$, and differential output state + * @f$y@f$, the DAE residual, derivative, and analytic Jacobian use + * @f[ + * T_1\dot y=K(u+b)-y. + * @f] + * Thus it is the existing GridDyn first-order lag/transducer primitive; it + * does not include a lead numerator time constant. At an equilibrium, + * @f$y=K(u+b)@f$, so desired-output initialization back-solves this relation. + * + * Parameters `t1` and `t` set @f$T_1@f$; inherited `k`/`gain`, `bias`, and + * output-limit parameters retain their GridBlock meanings. A time constant + * below the GridDyn numerical-resolution threshold selects the historical + * simplified gain mode. The solver path continues to use the equation above; + * the separate @ref step path uses the legacy local integration routine and is + * therefore not an exact sampled-data discretization. + */ class DelayBlock: public GridBlock { public: protected: - model_parameter mT1 = 0.1; //!< the time constant + model_parameter mT1 = 0.1; //!< Lag denominator time constant @f$T_1@f$. public: - //!< default constructor + /** @brief Construct a unity-gain lag with @f$T_1=0.1@f$. */ explicit DelayBlock(const std::string& objName = "delayBlock_#"); - /** alternate constructor to add in the time constant -@param[in] timeConstant the time constant -@param[in] objName the name of the block -*/ + /** + * @brief Construct a unity-gain lag. + * @param[in] timeConstant Lag time constant @f$T_1@f$. + * @param[in] objName Name of the block. + */ DelayBlock(double timeConstant, const std::string& objName = "delayBlock_#"); - /** alternate constructor to add in the time constant -@param[in] timeConstant the time constant -@param[in] gainValue the block gain -@param[in] objName the name of the object -*/ + /** + * @brief Construct a lag with explicit gain. + * @param[in] timeConstant Lag time constant @f$T_1@f$. + * @param[in] gainValue Steady-state gain @f$K@f$. + * @param[in] objName Name of the block. + */ DelayBlock(double timeConstant, double gainValue, const std::string& objName = "delayBlock_#"); virtual CoreObject* clone(CoreObject* obj = nullptr) const override; diff --git a/src/griddyn/blocks/DerivativeBlock.h b/src/griddyn/blocks/DerivativeBlock.h index 96c1b0b98..de479b05e 100644 --- a/src/griddyn/blocks/DerivativeBlock.h +++ b/src/griddyn/blocks/DerivativeBlock.h @@ -10,10 +10,16 @@ #include namespace griddyn::blocks { -/** @brief class implementing a derivative -block implementing \f$H(S)=\frac{K s}{1+T_1 s}\f$ -if the time constant is very small it reverts to the basic block -*/ +/** + * @brief First-order filtered differentiator. + * + * The block realizes @f$H(s)=Ks/(1+T_1s)@f$ through a filtered input state + * @f$z@f$: @f$T_1\dot z=K(u+b)-z@f$ and exposed output @f$y=\dot z@f$. + * Thus its DC output is zero and a constant input initializes the derivative + * to zero. `t1` or `t` sets @f$T_1@f$; inherited `k`/`gain`, `bias`, and + * output limits retain their GridBlock meanings. The solver uses the + * analytic residual/Jacobian; local stepping uses the legacy integration path. + */ class DerivativeBlock: public GridBlock { protected: model_parameter mT1 = 0.1; //!< delay time constant for the derivative filtering operation diff --git a/src/griddyn/blocks/FilteredDerivativeBlock.h b/src/griddyn/blocks/FilteredDerivativeBlock.h index 6867135a6..3e742cf73 100644 --- a/src/griddyn/blocks/FilteredDerivativeBlock.h +++ b/src/griddyn/blocks/FilteredDerivativeBlock.h @@ -10,10 +10,16 @@ #include namespace griddyn::blocks { -/** @brief class implementing a derivative -block implementing \f$H(S)=\frac{K s}{1+T_1 s} \frac{1}{1+T_2 s}\f$ - -*/ +/** + * @brief Two-stage filtered differentiator. + * + * This block realizes @f$H(s)=Ks/[(1+T_1s)(1+T_2s)]@f$. Its first state is + * the output of @f$T_1\dot z=K(u+b)-z@f$; the exposed output follows + * @f$T_2\dot y=\dot z-y@f$. Consequently a constant input has zero + * steady-state output. `t1` sets the pre-derivative filter and `t2` the + * output filter. Inherited output/ramp limits act on the exposed differential + * output and use normal GridBlock root handling. + */ class FilteredDerivativeBlock: public GridBlock { protected: model_parameter mT1 = 0.1; //!< delay time constant for the derivative filtering operation diff --git a/src/griddyn/blocks/FunctionBlock.cpp b/src/griddyn/blocks/FunctionBlock.cpp index 27507efe6..84966fc71 100644 --- a/src/griddyn/blocks/FunctionBlock.cpp +++ b/src/griddyn/blocks/FunctionBlock.cpp @@ -6,6 +6,7 @@ #include "FunctionBlock.h" +#include "core/CoreExceptions.h" #include "core/CoreObjectTemplates.hpp" #include "gmlc/utilities/stringOps.h" #include "gmlc/utilities/vectorOps.hpp" @@ -161,9 +162,7 @@ void FunctionBlock::setFunction(const std::string& functionName) mBinaryFunctionPtr = binaryFunctionPtr; opFlags.set(USES_CONSTANT_ARG); } else { - mFunctionPtr = nullptr; - mDerivativeFunctionPtr = nullptr; - mBinaryFunctionPtr = nullptr; + throw InvalidParameterValue("unknown function block function: " + functionName); } } diff --git a/src/griddyn/blocks/FunctionBlock.h b/src/griddyn/blocks/FunctionBlock.h index e4f1780dc..71fd65ea9 100644 --- a/src/griddyn/blocks/FunctionBlock.h +++ b/src/griddyn/blocks/FunctionBlock.h @@ -11,9 +11,17 @@ #include namespace griddyn::blocks { -/** @brief class implementing a function operation on the input -a wide assortment of functions are available including trig, logs, and other common math -operations*/ +/** + * @brief Algebraic wrapper around a supported unary or binary math function. + * + * Unary functions evaluate @f$y=K f(G(u+b))@f$. Binary functions evaluate + * @f$y=K f(G(u+b),a)@f$, where `arg` supplies the constant second argument. + * `function` and `func` select the function; an unknown name throws + * InvalidParameterValue rather than leaving a null callable. The analytic + * Jacobian is available for supported unary functions. Binary functions use + * the function interpreter's two-argument behavior and should be selected + * only when its derivative contract is suitable for the containing model. + */ class FunctionBlock: public GridBlock { public: //!< flags for function block diff --git a/src/griddyn/blocks/IntegralBlock.h b/src/griddyn/blocks/IntegralBlock.h index ebb8fa03d..ce391c168 100644 --- a/src/griddyn/blocks/IntegralBlock.h +++ b/src/griddyn/blocks/IntegralBlock.h @@ -10,9 +10,16 @@ #include namespace griddyn::blocks { -/** @brief class implementing an integral block -computes the integral of the input -*/ +/** + * @brief Differential integrator with optional GridBlock limits. + * + * For output state @f$y@f$, the governing equation is + * @f$\dot y=K(u+b)@f$. `iv`/`initial_value` specifies the state used when + * no desired output is supplied; `t` is accepted as an inverse gain, setting + * @f$K=1/t@f$. Inherited value and ramp limits use GridBlock's root and + * anti-windup behavior. The local stepping path uses trapezoidal input + * integration. + */ class IntegralBlock: public GridBlock { public: protected: diff --git a/src/griddyn/blocks/LeadLag.h b/src/griddyn/blocks/LeadLag.h new file mode 100644 index 000000000..55249f618 --- /dev/null +++ b/src/griddyn/blocks/LeadLag.h @@ -0,0 +1,124 @@ +/* + * Copyright (c) 2014-2026, Lawrence Livermore National Security + * See the top-level NOTICE for additional details. All rights reserved. + * SPDX-License-Identifier: BSD-3-Clause + */ + +#pragma once + +/** + * @file LeadLag.h + * @brief Stateless first-order lead--lag equations reusable by control models. + */ + +#include + +namespace griddyn::blocks { +/** + * @brief Stateless equations for a first-order lead--lag section. + * + * The kernel realizes + * @f[ + * G(s) = K\frac{1 + T_a s}{1 + T_b s} + * @f] + * with the lag state @f$x@f$: + * @f[ + * T_b\dot{x}=u-x,\qquad + * y=K\left[x+\frac{T_a}{T_b}(u-x)\right]. + * @f] + * Keeping these equations independent of @ref GridBlock lets composite + * controller models reuse the same realization while retaining ownership of + * their aggregate solver states, roots, and limiter policy. + * + * Some controller specifications instead expose the output itself as the + * differential state. For that state @f$y@f$, the equivalent realization is + * @f[ T_b\dot{y}=u-y+T_a\dot{u}. @f] + * @ref outputStateDerivative provides that form for models, such as TGOV1, + * whose published equations use it directly. + */ +class LeadLagKernel { + public: + /** Construct a unity-gain lag section with @f$T_a=0@f$ and @f$T_b=1@f$. */ + constexpr LeadLagKernel() = default; + /** Construct the section with numerator time @p leadTime, denominator time + * @p lagTime, and gain @p gainValue. */ + constexpr LeadLagKernel(double leadTime, double lagTime, double gainValue = 1.0): + Ta(leadTime), Tb(lagTime), K(gainValue) + { + } + + /** + * @brief Set all transfer-function parameters. + * + * @param[in] leadTime Numerator time constant @f$T_a@f$; it may be + * negative when a published output-state model requires it. + * @param[in] lagTime Positive denominator time constant @f$T_b@f$. + * @param[in] gainValue Transfer-function gain @f$K@f$. + */ + constexpr void setParameters(double leadTime, double lagTime, double gainValue = 1.0) + { + Ta = leadTime; + Tb = lagTime; + K = gainValue; + } + /** @return true if all parameters are finite and @f$T_b>0@f$. */ + [[nodiscard]] bool isValid() const + { + return std::isfinite(Ta) && std::isfinite(Tb) && std::isfinite(K) && (Tb > 0.0); + } + /** @return the lead numerator time constant. */ + [[nodiscard]] constexpr double leadTime() const { return Ta; } + /** @return the lag denominator time constant. */ + [[nodiscard]] constexpr double lagTime() const { return Tb; } + /** @return the section gain. */ + [[nodiscard]] constexpr double gain() const { return K; } + + /** + * @brief Evaluate the unbounded lag-state output @f$y@f$. + * @param[in] input Input @f$u@f$. + * @param[in] state Lag state @f$x@f$. + */ + [[nodiscard]] double output(double input, double state) const + { + return K * (state + ((Ta / Tb) * (input - state))); + } + /** + * @brief Evaluate the lag-state derivative @f$\dot{x}@f$. + * @param[in] input Input @f$u@f$. + * @param[in] state Lag state @f$x@f$. + */ + [[nodiscard]] double derivative(double input, double state) const + { + return (input - state) / Tb; + } + /** + * @brief Evaluate @f$\dot y=(u-y+T_a\dot u)/T_b@f$. + * + * This is the equivalent output-state realization, used when a model + * explicitly owns @f$y@f$ rather than the lag state @f$x@f$. + * @param[in] input Input @f$u@f$. + * @param[in] outputState Output state @f$y@f$. + * @param[in] inputDerivative Input derivative @f$\dot u@f$. + */ + [[nodiscard]] + double outputStateDerivative(double input, double outputState, double inputDerivative) const + { + return (input - outputState + (Ta * inputDerivative)) / Tb; + } + /** @return @f$\partial y/\partial u@f$. */ + [[nodiscard]] double outputInputJacobian() const { return K * Ta / Tb; } + /** @return @f$\partial y/\partial x@f$. */ + [[nodiscard]] double outputStateJacobian() const { return K * (1.0 - (Ta / Tb)); } + /** @return @f$\partial\dot{x}/\partial u@f$. */ + [[nodiscard]] double derivativeInputJacobian() const { return 1.0 / Tb; } + /** @return @f$\partial\dot{x}/\partial x@f$. */ + [[nodiscard]] double derivativeStateJacobian() const { return -1.0 / Tb; } + /** @return the coefficient of @f$\dot u@f$ in the output-state realization. */ + [[nodiscard]] double outputStateInputDerivativeJacobian() const { return Ta / Tb; } + + private: + double Ta = 0.0; + double Tb = 1.0; + double K = 1.0; +}; +} // namespace griddyn::blocks diff --git a/src/griddyn/blocks/LeadLagBlock.cpp b/src/griddyn/blocks/LeadLagBlock.cpp new file mode 100644 index 000000000..bfad0ce09 --- /dev/null +++ b/src/griddyn/blocks/LeadLagBlock.cpp @@ -0,0 +1,206 @@ +/* + * Copyright (c) 2014-2026, Lawrence Livermore National Security + * See the top-level NOTICE for additional details. All rights reserved. + * SPDX-License-Identifier: BSD-3-Clause + */ + +#include "LeadLagBlock.h" + +#include "ValueLimiter.h" +#include "core/CoreExceptions.h" +#include "core/CoreObjectTemplates.hpp" +#include "utilities/MatrixData.hpp" +#include + +namespace griddyn::blocks { +LeadLagBlock::LeadLagBlock(const std::string& objName): GridBlock(objName) +{ + opFlags.set(USE_STATE); +} + +LeadLagBlock::LeadLagBlock(double lagTime, double leadTime, const std::string& objName): + GridBlock(objName), section(leadTime, lagTime) +{ + opFlags.set(USE_STATE); +} + +LeadLagBlock::LeadLagBlock(double lagTime, + double leadTime, + double gainValue, + const std::string& objName): + GridBlock(gainValue, objName), section(leadTime, lagTime, gainValue) +{ + opFlags.set(USE_STATE); +} + +CoreObject* LeadLagBlock::clone(CoreObject* obj) const +{ + auto* clone = cloneBase(this, obj); + if (clone != nullptr) { + clone->section = section; + } + return (clone == nullptr) ? obj : clone; +} + +void LeadLagBlock::validateParameters() const +{ + if (!section.isValid()) { + throw InvalidParameterValue("lead-lag gain or time constants"); + } +} + +void LeadLagBlock::dynObjectInitializeA(CoreTime time0, std::uint32_t flags) +{ + validateParameters(); + GridBlock::dynObjectInitializeA(time0, flags); + // GridBlock supplies the algebraic output. This block adds the one + // differential lag state used by the reusable lead-lag realization. + ++offsets.local().local.diffSize; + offsets.local().local.jacSize += 4; +} + +void LeadLagBlock::dynObjectInitializeB(const IOdata& inputs, + const IOdata& desiredOutput, + IOdata& fieldSet) +{ + fieldSet.resize(1); + const index_t lagState = offsets.local().local.algSize; + if (desiredOutput.empty()) { + const double input = inputs[0] + bias; + m_state[lagState] = input; + m_state[limiter_alg] = section.output(input, input); + if (opFlags[USE_BLOCK_LIMITS]) { + m_state[0] = vLimiter->clampOutput(m_state[limiter_alg]); + } + fieldSet[0] = m_state[0]; + prevInput = input; + return; + } + + const double output = + opFlags[USE_BLOCK_LIMITS] ? vLimiter->clampOutput(desiredOutput[0]) : desiredOutput[0]; + const double input = output / section.gain(); + m_state[lagState] = input; + m_state[limiter_alg] = output; + if (opFlags[USE_BLOCK_LIMITS]) { + m_state[0] = output; + } + fieldSet[0] = input - bias; + prevInput = input; +} + +void LeadLagBlock::blockDerivative(double input, + double /*didt*/, + const StateData& stateDataValue, + double deriv[], + const SolverMode& solverModeValue) +{ + const auto locations = offsets.getLocations(stateDataValue, deriv, solverModeValue, this); + locations.destDiffLoc[0] = section.derivative(input + bias, locations.diffStateLoc[0]); +} + +void LeadLagBlock::blockAlgebraicUpdate(double input, + const StateData& stateDataValue, + double update[], + const SolverMode& solverModeValue) +{ + const auto locations = offsets.getLocations(stateDataValue, update, solverModeValue, this); + locations.destLoc[limiter_alg] = section.output(input + bias, locations.diffStateLoc[0]); + if (limiter_alg > 0) { + GridBlock::blockAlgebraicUpdate(input, stateDataValue, update, solverModeValue); + } +} + +void LeadLagBlock::blockJacobianElements(double /*input*/, + double /*didt*/, + const StateData& stateDataValue, + MatrixData& matrixDataValue, + index_t argLoc, + const SolverMode& solverModeValue) +{ + const auto locations = offsets.getLocations(stateDataValue, solverModeValue, this); + if (hasAlgebraic(solverModeValue)) { + const auto outputLocation = locations.algOffset + limiter_alg; + matrixDataValue.assign(outputLocation, outputLocation, -1.0); + matrixDataValue.assignCheckCol(outputLocation, argLoc, section.outputInputJacobian()); + if (hasDifferential(solverModeValue)) { + matrixDataValue.assign(outputLocation, + locations.diffOffset, + section.outputStateJacobian()); + } + } + if (hasDifferential(solverModeValue)) { + matrixDataValue.assignCheckCol(locations.diffOffset, + argLoc, + section.derivativeInputJacobian()); + matrixDataValue.assign(locations.diffOffset, + locations.diffOffset, + section.derivativeStateJacobian() - stateDataValue.cj); + } + if ((limiter_alg > 0) && hasAlgebraic(solverModeValue)) { + GridBlock::blockJacobianElements( + 0.0, 0.0, stateDataValue, matrixDataValue, argLoc, solverModeValue); + } +} + +double LeadLagBlock::step(CoreTime time, double inputValue) +{ + const double input = inputValue + bias; + const double timeStep = time - prevTime; + const index_t lagState = offsets.local().local.algSize; + if (timeStep > 0.0) { + // Exact zero-order-hold solution of T_b xdot=u-x. Solver-based + // execution uses the residual above; this path serves local stepping. + const double decay = std::exp(-timeStep / section.lagTime()); + m_state[lagState] = input + ((m_state[lagState] - input) * decay); + } + m_state[limiter_alg] = section.output(input, m_state[lagState]); + prevInput = input; + if (opFlags[USE_BLOCK_LIMITS]) { + return GridBlock::step(time, inputValue); + } + prevTime = time; + m_output = m_state[0]; + return m_output; +} + +void LeadLagBlock::set(std::string_view param, std::string_view val) +{ + GridBlock::set(param, val); +} + +void LeadLagBlock::set(std::string_view param, double val, units::unit unitType) +{ + if (param == "ta") { + if (!std::isfinite(val)) { + throw InvalidParameterValue("lead-lag Ta must be finite"); + } + section.setParameters(val, section.lagTime(), section.gain()); + } else if ((param == "tb") || (param == "t")) { + if (!std::isfinite(val) || (val <= 0.0)) { + throw InvalidParameterValue("lead-lag Tb must be positive and finite"); + } + section.setParameters(section.leadTime(), val, section.gain()); + } else if ((param == "k") || (param == "gain")) { + if (!std::isfinite(val)) { + throw InvalidParameterValue("lead-lag gain must be finite"); + } + GridBlock::set(param, val, unitType); + section.setParameters(section.leadTime(), section.lagTime(), val); + } else { + GridBlock::set(param, val, unitType); + } +} + +stringVec LeadLagBlock::localStateNames() const +{ + stringVec names(stateSize(cLocalSolverMode)); + index_t index = 0; + if (opFlags[USE_BLOCK_LIMITS]) { + names[index++] = "output"; + } + names[index++] = "unlimited_output"; + names[index] = "lag_state"; + return names; +} +} // namespace griddyn::blocks diff --git a/src/griddyn/blocks/LeadLagBlock.h b/src/griddyn/blocks/LeadLagBlock.h new file mode 100644 index 000000000..f961bd24e --- /dev/null +++ b/src/griddyn/blocks/LeadLagBlock.h @@ -0,0 +1,94 @@ +/* + * Copyright (c) 2014-2026, Lawrence Livermore National Security + * See the top-level NOTICE for additional details. All rights reserved. + * SPDX-License-Identifier: BSD-3-Clause + */ + +#pragma once + +/** + * @file LeadLagBlock.h + * @brief GridBlock implementation of a first-order lead--lag transfer function. + */ + +#include "../Block.h" +#include "LeadLag.h" +#include + +namespace griddyn::blocks { +/** + * @brief GridBlock wrapper for a first-order lead--lag section. + * + * @f[ + * G(s)=K\frac{1+T_a s}{1+T_b s},\qquad + * T_b\dot{x}=u+b-x,\qquad + * y=K\left[x+\frac{T_a}{T_b}(u+b-x)\right]. + * @f] + * + * The output is algebraic and the lag state is differential. Configured + * GridBlock output limits apply to the exposed output only; they do not alter + * the linear section's state equation. A non-positive denominator time is + * invalid rather than silently changing this block into a gain, because that + * would change its state and solver contract after construction. + * + * The block may be created through the `block` factory as `leadlag` or + * `lead_lag`. Parameters `ta` and `tb` set @f$T_a@f$ and @f$T_b@f$; + * `t` aliases `tb`; inherited `k`/`gain` and `bias` set @f$K@f$ and @f$b@f$. + * At equilibrium, @f$x=u+b@f$ and @f$y=K(u+b)@f$. Initialization from a + * desired output back-solves that steady-state relation. The local + * @ref step method uses the exact zero-order-hold solution for @f$x@f$; + * solver execution uses the residual and analytic Jacobian above. + */ +class LeadLagBlock: public GridBlock { + public: + /** @brief Construct a unity-gain block with @f$T_a=0@f$, @f$T_b=1@f$. */ + explicit LeadLagBlock(const std::string& objName = "leadLagBlock_#"); + /** + * @brief Construct a unity-gain lead--lag block. + * @param[in] lagTime Denominator time constant @f$T_b@f$. + * @param[in] leadTime Numerator time constant @f$T_a@f$. + * @param[in] objName Name of the block. + */ + LeadLagBlock(double lagTime, double leadTime, const std::string& objName = "leadLagBlock_#"); + /** + * @brief Construct a lead--lag block with explicit gain. + * @param[in] lagTime Denominator time constant @f$T_b@f$. + * @param[in] leadTime Numerator time constant @f$T_a@f$. + * @param[in] gainValue Transfer-function gain @f$K@f$. + * @param[in] objName Name of the block. + */ + LeadLagBlock(double lagTime, + double leadTime, + double gainValue, + const std::string& objName = "leadLagBlock_#"); + CoreObject* clone(CoreObject* obj = nullptr) const override; + + void dynObjectInitializeA(CoreTime time0, std::uint32_t flags) override; + void dynObjectInitializeB(const IOdata& inputs, + const IOdata& desiredOutput, + IOdata& fieldSet) override; + void set(std::string_view param, std::string_view val) override; + void set(std::string_view param, double val, units::unit unitType = units::defunit) override; + void blockDerivative(double input, + double didt, + const StateData& stateDataValue, + double deriv[], + const SolverMode& solverModeValue) override; + void blockAlgebraicUpdate(double input, + const StateData& stateDataValue, + double update[], + const SolverMode& solverModeValue) override; + void blockJacobianElements(double input, + double didt, + const StateData& stateDataValue, + MatrixData& matrixDataValue, + index_t argLoc, + const SolverMode& solverModeValue) override; + double step(CoreTime time, double input) override; + stringVec localStateNames() const override; + + private: + LeadLagKernel section; + void validateParameters() const; +}; +} // namespace griddyn::blocks diff --git a/src/griddyn/blocks/LutBlock.cpp b/src/griddyn/blocks/LutBlock.cpp index 45b1f26f2..27835a5d5 100644 --- a/src/griddyn/blocks/LutBlock.cpp +++ b/src/griddyn/blocks/LutBlock.cpp @@ -6,12 +6,14 @@ #include "LutBlock.h" +#include "ValueLimiter.h" +#include "core/CoreExceptions.h" #include "core/CoreObjectTemplates.hpp" #include "gmlc/utilities/TimeSeries.hpp" #include "gmlc/utilities/stringConversion.h" -#include "gmlc/utilities/vectorOps.hpp" #include "utilities/MatrixData.hpp" #include +#include #include #include @@ -27,11 +29,6 @@ CoreObject* LutBlock::clone(CoreObject* obj) const return obj; } nobj->lut = lut; - nobj->b = b; - nobj->m = m; - nobj->vlower = vlower; - nobj->vupper = vupper; - nobj->lindex = lindex; return nobj; } @@ -40,13 +37,45 @@ void LutBlock::dynObjectInitializeB(const IOdata& inputs, const IOdata& desiredOutput, IOdata& fieldSet) { + validateTable(lut); + fieldSet.resize(1); + double input = inputs.empty() ? 0.0 : inputs[0] + bias; + if (desiredOutput.empty()) { - m_state[limiter_alg] = K * computeValue(inputs[0] + bias); - GridBlock::dynObjectInitializeB(inputs, desiredOutput, fieldSet); + if (inputs.empty()) { + throw InvalidParameterValue("LUT initialization requires an input or desired output"); + } } else { - // TODO(pt): figure out how to invert the lookup table - GridBlock::dynObjectInitializeB(inputs, desiredOutput, fieldSet); + if (!std::isfinite(K) || (std::abs(K) < kMin_Res)) { + throw InvalidParameterValue( + "LUT desired-output initialization requires a finite nonzero gain"); + } + const double limitedOutput = opFlags[USE_BLOCK_LIMITS] ? + std::clamp(desiredOutput[0], static_cast(Omin), static_cast(Omax)) : + desiredOutput[0]; + input = inverseValue(limitedOutput / K); + fieldSet[0] = input - bias; + } + + const double output = K * evaluate(input).value; + if (opFlags[USE_BLOCK_LIMITS]) { + m_state[limiter_alg] = output; + GridBlock::rootCheck({input - bias}, + emptyStateData, + cLocalSolverMode, + CheckLevel::REVERSABLE_ONLY); + m_state[0] = vLimiter->clampOutput(m_state[limiter_alg]); + } else { + m_state[0] = output; + m_output = output; + } + if (desiredOutput.empty()) { + fieldSet[0] = opFlags[USE_BLOCK_LIMITS] ? m_state[0] : output; + } + if (opFlags[USE_BLOCK_LIMITS]) { + m_output = m_state[0]; } + prevInput = input; } void LutBlock::blockAlgebraicUpdate(double input, @@ -55,7 +84,7 @@ void LutBlock::blockAlgebraicUpdate(double input, const SolverMode& sMode) { auto offset = offsets.getAlgOffset(sMode) + limiter_alg; - update[offset] = K * computeValue(input + bias); + update[offset] = K * evaluate(input + bias).value; if (limiter_alg > 0) { GridBlock::blockAlgebraicUpdate(input, stateDataValue, update, sMode); return; @@ -72,7 +101,7 @@ void LutBlock::blockJacobianElements(double input, auto offset = offsets.getAlgOffset(sMode) + limiter_alg; // use the md.assign Macro defined in basicDefs // md.assign(arrayIndex, RowIndex, ColIndex, value) - matrixDataValue.assignCheckCol(offset, argLoc, K * m); + matrixDataValue.assignCheckCol(offset, argLoc, K * evaluate(input + bias).slope); matrixDataValue.assign(offset, offset, -1); if (limiter_alg > 0) { GridBlock::blockJacobianElements( @@ -86,36 +115,42 @@ void LutBlock::set(std::string_view param, std::string_view val) using gmlc::utilities::str2vector; if (param == "lut") { const auto vectorData = str2vector(std::string{val}, -kBigNum, ";,:"); - lut.clear(); - lut.emplace_back(-kBigNum, 0.0); - lut.emplace_back(kBigNum, 0.0); - for (size_t mm = 0; mm < vectorData.size(); mm += 2) { - lut.emplace_back(vectorData[mm], vectorData[mm + 1]); + if ((vectorData.size() % 2) != 0) { + throw InvalidParameterValue("LUT table requires x,y point pairs"); + } + std::vector> table; + table.reserve(vectorData.size() / 2); + for (size_t index = 0; index < vectorData.size(); index += 2) { + table.emplace_back(vectorData[index], vectorData[index + 1]); } - std::sort(lut.begin(), lut.end()); - lut[0].second = lut[1].second; - (*lut.end()).second = (*(lut.end() - 1)).second; + std::sort(table.begin(), table.end()); + validateTable(table); + lut = std::move(table); } else if (param == "element") { const auto vectorData = str2vector(std::string{val}, -kBigNum, ";,:"); - for (size_t mm = 0; mm < vectorData.size(); mm += 2) { - lut.emplace_back(vectorData[mm], vectorData[mm + 1]); + if ((vectorData.size() % 2) != 0) { + throw InvalidParameterValue("LUT table requires x,y point pairs"); + } + auto table = lut; + table.reserve(table.size() + (vectorData.size() / 2)); + for (size_t index = 0; index < vectorData.size(); index += 2) { + table.emplace_back(vectorData[index], vectorData[index + 1]); } - std::sort(lut.begin(), lut.end()); - lut[0].second = lut[1].second; - (*lut.end()).second = (*(lut.end() - 1)).second; + std::sort(table.begin(), table.end()); + validateTable(table); + lut = std::move(table); } else if (param == "file") { const gmlc::utilities::TimeSeries timeSeries(std::string{val}); - lut.clear(); - lut.emplace_back(-kBigNum, 0.0); - lut.emplace_back(kBigNum, 0.0); + std::vector> table; + table.reserve(timeSeries.size()); for (gmlc::utilities::fsize_t pointIndex = 0; pointIndex < timeSeries.size(); ++pointIndex) { - lut.emplace_back(timeSeries.time(pointIndex), timeSeries.data(pointIndex)); + table.emplace_back(timeSeries.time(pointIndex), timeSeries.data(pointIndex)); } - std::sort(lut.begin(), lut.end()); - lut[0].second = lut[1].second; - (*lut.end()).second = (*(lut.end() - 1)).second; + std::sort(table.begin(), table.end()); + validateTable(table); + lut = std::move(table); } else { GridBlock::set(param, val); } @@ -123,49 +158,101 @@ void LutBlock::set(std::string_view param, std::string_view val) void LutBlock::set(std::string_view param, double val, units::unit unitType) { - if (param.empty() || param[0] == '#') { - } else { + if (!param.empty() && param[0] != '#') { GridBlock::set(param, val, unitType); } } double LutBlock::step(CoreTime time, double input) { - m_state[limiter_alg] = K * computeValue(input + bias); + const double output = K * evaluate(input + bias).value; if (limiter_alg > 0) { + m_state[limiter_alg] = output; GridBlock::step(time, input); } else { - m_output = m_state[0]; + m_state[0] = output; + m_output = output; prevTime = time; } - return m_state[0]; + return (limiter_alg > 0) ? m_state[0] : output; } -double LutBlock::computeValue(double input) +void LutBlock::validateTable(const std::vector>& table) { - if (input > vupper) { - ++lindex; - auto lower = std::lower_bound(lut.begin() + lindex, lut.end(), std::make_pair(input, 0.0)); - auto upper = lower; - ++upper; - lindex = static_cast(upper - lut.begin()); - vlower = lower->first; - vupper = upper->first; - m = (upper->second - lower->second) / (vupper - vlower); - b = lower->second; - } else if (input < vlower) { - --lindex; - while (lut[lindex].first > input) { - --lindex; + if (table.empty()) { + throw InvalidParameterValue("LUT table requires at least one point"); + } + for (size_t index = 0; index < table.size(); ++index) { + if (!std::isfinite(table[index].first) || !std::isfinite(table[index].second)) { + throw InvalidParameterValue("LUT table values must be finite"); + } + if ((index > 0) && !(table[index].first > table[index - 1].first)) { + throw InvalidParameterValue("LUT table abscissas must be strictly increasing"); } - vlower = lut[lindex - 1].first; - vupper = lut[lindex].first; - m = (lut[lindex].second - lut[lindex - 1].second) / (vupper - vlower); - b = lut[lindex - 1].second; } - return (((input - vlower) * m) + b); +} + +LutBlock::LookupResult LutBlock::evaluate(double input) const +{ + validateTable(lut); + if (lut.size() == 1 || input <= lut.front().first) { + return {lut.front().second, 0.0}; + } + if (input >= lut.back().first) { + return {lut.back().second, 0.0}; + } + + const auto upper = + std::upper_bound(lut.begin(), lut.end(), input, [](double value, const auto& point) { + return value < point.first; + }); + const auto lower = std::prev(upper); + const double slope = (upper->second - lower->second) / (upper->first - lower->first); + return {lower->second + ((input - lower->first) * slope), slope}; +} + +double LutBlock::inverseValue(double value) const +{ + validateTable(lut); + bool nondecreasing = true; + bool nonincreasing = true; + for (size_t index = 1; index < lut.size(); ++index) { + nondecreasing = nondecreasing && (lut[index].second >= lut[index - 1].second); + nonincreasing = nonincreasing && (lut[index].second <= lut[index - 1].second); + } + if (!nondecreasing && !nonincreasing) { + throw InvalidParameterValue("LUT desired-output initialization requires a monotonic table"); + } + + const double minimum = std::min(lut.front().second, lut.back().second); + const double maximum = std::max(lut.front().second, lut.back().second); + if ((value < minimum) || (value > maximum)) { + throw InvalidParameterValue("LUT desired output is outside the table range"); + } + if (lut.size() == 1) { + return lut.front().first; + } + for (size_t index = 1; index < lut.size(); ++index) { + const auto& lower = lut[index - 1]; + const auto& upper = lut[index]; + if ((value >= std::min(lower.second, upper.second)) && + (value <= std::max(lower.second, upper.second))) { + if (upper.second == lower.second) { + return lower.first; + } + return lower.first + + ((value - lower.second) * (upper.first - lower.first) / + (upper.second - lower.second)); + } + } + return lut.back().first; +} + +double LutBlock::computeValue(double input) const +{ + return evaluate(input).value; } } // namespace griddyn::blocks diff --git a/src/griddyn/blocks/LutBlock.h b/src/griddyn/blocks/LutBlock.h index fff163d6d..e2c9c8f5e 100644 --- a/src/griddyn/blocks/LutBlock.h +++ b/src/griddyn/blocks/LutBlock.h @@ -8,20 +8,45 @@ #include "../Block.h" #include +#include #include namespace griddyn::blocks { -/** @brief lookup table block*/ +/** + * @brief Static piecewise-linear lookup-table block. + * + * Given strictly increasing abscissas @f$x_i@f$ and tabulated values @f$y_i@f$, + * this block evaluates + * @f[ + * y = K f(u + b), + * @f] + * where @f$f@f$ is linear on each interval @f$[x_i,x_{i+1}]@f$ and is held + * at the nearest endpoint outside the table domain. Its analytic input + * Jacobian is @f$K (y_{i+1}-y_i)/(x_{i+1}-x_i)@f$ in an interval and zero + * in either endpoint-held region. + * + * The table may be set with @c lut (replacing all points), @c element + * (appending points), or @c file. A desired-output initialization is + * available only for a monotonic table with a finite nonzero gain, since a + * non-monotonic table does not have an unambiguous inverse. + */ class LutBlock: public GridBlock { public: private: std::vector> lut; //!< the lookup table - double b = 0; //!< the intercept of the interpolation function of the current lookup section - double m = 0; //!< the slope of the interpolation function of the current lookup section - double vlower = -kBigNum; //!< the lower value of the current lookup table section - double vupper = kBigNum; //!< the upper value of the current lookup table section - int lindex = -1; //!< the index of the current lookup table section - // NOTE: extra 4 bytes here + + struct LookupResult { + double value; + double slope; + }; + + /** @brief Check that a table can be evaluated without ambiguity or division by zero. */ + static void validateTable(const std::vector>& table); + /** @brief Evaluate the table and its local derivative without mutable solver-side caching. */ + [[nodiscard]] LookupResult evaluate(double input) const; + /** @brief Return the input that produces a table value for monotonic-table initialization. */ + [[nodiscard]] double inverseValue(double value) const; + public: explicit LutBlock(const std::string& objName = "lutBlock_#"); virtual CoreObject* clone(CoreObject* obj = nullptr) const override; @@ -48,6 +73,7 @@ class LutBlock: public GridBlock { const SolverMode& sMode) override; virtual double step(CoreTime time, double input) override; // virtual void setTime(CoreTime time){prevTime=time;}; - double computeValue(double input); + /** @brief Return the endpoint-clamped, piecewise-linear table value. */ + [[nodiscard]] double computeValue(double input) const; }; } // namespace griddyn::blocks diff --git a/src/griddyn/blocks/NullBlock.h b/src/griddyn/blocks/NullBlock.h index 9dd0773a2..e59de2ad8 100644 --- a/src/griddyn/blocks/NullBlock.h +++ b/src/griddyn/blocks/NullBlock.h @@ -11,7 +11,13 @@ #include namespace griddyn::blocks { -/** @brief class defining a null block meaning input==output +/** + * @brief Identity/pass-through block with no solver states or roots. + * + * The block implements @f$y=u@f$ and, when configured with a differential + * input, passes its input derivative through as well. It deliberately has no + * GridBlock gain, bias, limiter, power-flow, or dynamic operations. It is a + * programmatic composition helper and is intentionally not factory-registered. */ class NullBlock final: public GridBlock { public: diff --git a/src/griddyn/blocks/PidBlock.h b/src/griddyn/blocks/PidBlock.h index f9c6a9c36..c4abea9cc 100644 --- a/src/griddyn/blocks/PidBlock.h +++ b/src/griddyn/blocks/PidBlock.h @@ -10,9 +10,17 @@ #include namespace griddyn::blocks { -/** @brief class implementing a PID controller -the derivative operator has a prefilter operation on it with a time constant T1 and the output has a -delay of Td*/ +/** + * @brief Filtered PID controller with a first-order output stage. + * + * With error @f$e=u+b@f$, integral state @f$i@f$, derivative-filter state + * @f$d@f$, and output state @f$y@f$, the model uses + * @f[\dot i=I e,\qquad T_1\dot d=D e-d,\qquad + * T_d\dot y=K(Pe+\dot d+i)-y.\f] + * The derivative path is omitted when `d` is zero. Parameters `p`, `i`, + * `d`, `t1`, `td`, and `iv`/`initial_value` set the displayed quantities; + * inherited GridBlock limits apply to the output state. + */ class PidBlock: public GridBlock { public: protected: diff --git a/src/griddyn/blocks/RampLimiter.h b/src/griddyn/blocks/RampLimiter.h index 24c1955c7..73da3eea7 100644 --- a/src/griddyn/blocks/RampLimiter.h +++ b/src/griddyn/blocks/RampLimiter.h @@ -8,12 +8,19 @@ #include namespace griddyn::blocks { -/** class that limits the rate of change of a value between an upper and lower limit and - * maintains state of whether it is clamping or not +/** + * @brief Stateful rate limiter used by dynamic control blocks. + * + * For requested rate @f$r@f$, the limiting law is + * @f$\hat r=\min(\max(r,r_{min}),r_{max})@f$. Unlike the stateless + * @ref clampOutputRamp helper, @ref output and @ref deriv retain engagement + * state. While engaged their input derivative is zero, providing the Jacobian + * behavior needed for root-based anti-windup. `resetLevel` supplies the + * hysteresis margin used before a held limiter releases. */ class RampLimiter { private: - double minRamp = std::numeric_limits::min(); //!< the minimum ramp + double minRamp = std::numeric_limits::lowest(); //!< the minimum ramp double maxRamp = std::numeric_limits::max(); //!< the maximum ramp to allow double resetLevel = 0.0; //!< the level by which an input ramp has to be below the limits to reset diff --git a/src/griddyn/blocks/TransferFunctionBlock.cpp b/src/griddyn/blocks/TransferFunctionBlock.cpp index 1c0b477d1..00363ce9b 100644 --- a/src/griddyn/blocks/TransferFunctionBlock.cpp +++ b/src/griddyn/blocks/TransferFunctionBlock.cpp @@ -6,288 +6,479 @@ #include "TransferFunctionBlock.h" +#include "ValueLimiter.h" #include "core/CoreExceptions.h" #include "core/CoreObjectTemplates.hpp" #include "gmlc/utilities/stringConversion.h" -#include "gmlc/utilities/vectorOps.hpp" +#include "gmlc/utilities/stringOps.h" #include "utilities/MatrixData.hpp" +#include +#include #include #include #include namespace griddyn::blocks { TransferFunctionBlock::TransferFunctionBlock(const std::string& objName): - GridBlock(objName), a(2, 1), b(2, 0) + GridBlock(objName), a{1.0, 1.0}, b{1.0, 0.0} { - b[0] = 1; opFlags.set(USE_STATE); } -TransferFunctionBlock::TransferFunctionBlock(int order): a(order + 1, 1), b(order + 1, 0) +TransferFunctionBlock::TransferFunctionBlock(int orderValue): + GridBlock("transferFunctionBlock_#"), a(static_cast(std::max(0, orderValue)) + 1, 0.0), + b(a.size(), 0.0) { - if (a.empty()) { - a.push_back(1.0); - b.assign(1, 1.0); - } else { - b[0] = 1; - } + a.front() = 1.0; + a.back() = 1.0; + b.front() = 1.0; opFlags.set(USE_STATE); } TransferFunctionBlock::TransferFunctionBlock(std::vector acoef): - a(std::move(acoef)), b(a.size(), 0) + GridBlock("transferFunctionBlock_#"), a(std::move(acoef)), b(a.size(), 0.0) { if (a.empty()) { a.push_back(1.0); b.assign(1, 1.0); } else { - b[0] = 1; + b.front() = 1.0; } opFlags.set(USE_STATE); } TransferFunctionBlock::TransferFunctionBlock(std::vector acoef, std::vector bcoef): - a(std::move(acoef)), b(std::move(bcoef)) + GridBlock("transferFunctionBlock_#"), a(std::move(acoef)), b(std::move(bcoef)) { if (a.empty()) { a.push_back(1.0); } - b.resize(a.size(), 0); + b.resize(a.size(), 0.0); opFlags.set(USE_STATE); } CoreObject* TransferFunctionBlock::clone(CoreObject* obj) const { - auto* nobj = cloneBase(this, obj); - if (nobj == nullptr) { + auto* clone = cloneBase(this, obj); + if (clone == nullptr) { return obj; } + clone->a = a; + clone->b = b; + return clone; +} - nobj->a = a; - nobj->b = b; - return nobj; +index_t TransferFunctionBlock::order() const +{ + return static_cast(a.size() - 1); } -// set up the number of states -void TransferFunctionBlock::dynObjectInitializeA(CoreTime time0, std::uint32_t flags) + +void TransferFunctionBlock::validateCoefficients() { - if (b.back() == 0) { - opFlags[DIFFERENTIAL_OUTPUT] = true; - extraOutputState = false; - } else { - extraOutputState = true; + if (a.empty()) { + throw InvalidParameterValue("transfer-function denominator"); + } + if (b.size() > a.size()) { + throw InvalidParameterValue("transfer-function numerator order"); + } + b.resize(a.size(), 0.0); + const bool invalidDenominator = + std::any_of(a.begin(), a.end(), [](double value) { return !std::isfinite(value); }) || + (std::abs(a.back()) < kMin_Res); + const bool invalidNumerator = + std::any_of(b.begin(), b.end(), [](double value) { return !std::isfinite(value); }); + if (invalidDenominator || invalidNumerator || !std::isfinite(K) || !std::isfinite(bias)) { + throw InvalidParameterValue("transfer-function coefficients, gain, or bias"); } +} + +void TransferFunctionBlock::dynObjectInitializeA(CoreTime time0, std::uint32_t flags) +{ + validateCoefficients(); + opFlags.reset(DIFFERENTIAL_OUTPUT); + // GridBlock's convenience use_limits flag also requests rate limits. This + // block has an algebraic output, so retain only its supported output limit. + opFlags.reset(USE_RAMP_LIMITS); GridBlock::dynObjectInitializeA(time0, flags); - offsets.local().local.jacSize += static_cast((3 * (a.size() - 2)) + 1); - offsets.local().local.diffSize += static_cast(a.size()) - 2; - if (extraOutputState) { - offsets.local().local.diffSize += 1; - offsets.local().local.jacSize += 3; + offsets.local().local.diffSize += order(); + offsets.local().local.jacSize += static_cast((3 * order()) + 3); +} + +double TransferFunctionBlock::rawOutput(double input, const double state[]) const +{ + const index_t stateCount = order(); + const double denominatorScale = a.back(); + if (stateCount == 0) { + return K * (b.front() / denominatorScale) * input; } + + const double directTerm = b.back() / denominatorScale; + double output = directTerm * input; + for (index_t index = 0; index < stateCount; ++index) { + output += ((b[index] / denominatorScale) - (directTerm * a[index] / denominatorScale)) * + state[index]; + } + return K * output; } -// initial conditions + +void TransferFunctionBlock::stateDerivative(double input, + const double state[], + double derivative[]) const +{ + const index_t stateCount = order(); + if (stateCount == 0) { + return; + } + for (index_t index = 0; index + 1 < stateCount; ++index) { + derivative[index] = state[index + 1]; + } + double finalDerivative = input; + const double denominatorScale = a.back(); + for (index_t index = 0; index < stateCount; ++index) { + finalDerivative -= (a[index] / denominatorScale) * state[index]; + } + derivative[stateCount - 1] = finalDerivative; +} + +double TransferFunctionBlock::externalStateOutput(double input, const double state[]) const +{ + return rawOutput(input, state); +} + +void TransferFunctionBlock::externalStateDerivative(double input, + const double state[], + double derivative[]) const +{ + stateDerivative(input, state, derivative); +} + void TransferFunctionBlock::dynObjectInitializeB(const IOdata& inputs, const IOdata& desiredOutput, IOdata& fieldSet) { - if (desiredOutput.empty()) { - // m_state[2] = (1.0 - m_T2 / m_T1) * (inputs[0] + bias); - m_state[1] = (inputs[0] + bias); - m_state[0] = m_state[1] * K; - if (opFlags[HAS_LIMITS]) { - GridBlock::rootCheck(inputs, - emptyStateData, - cLocalSolverMode, - CheckLevel::REVERSABLE_ONLY); - m_state[0] = gmlc::utilities::valLimit(m_state[0], Omin, Omax); + fieldSet.resize(1); + const index_t stateCount = order(); + const index_t rawOutputIndex = limiter_alg; + const index_t stateStart = limiter_alg + 1; + double input = inputs.empty() ? 0.0 : inputs[0] + bias; + + if (!desiredOutput.empty()) { + const double dcGain = K * b.front() / a.front(); + if (!std::isfinite(dcGain) || (std::abs(a.front()) < kMin_Res) || + (std::abs(dcGain) < kMin_Res)) { + throw InvalidParameterValue( + "transfer-function desired-output initialization requires finite nonzero DC gain"); } + const double limitedOutput = opFlags[USE_BLOCK_LIMITS] ? + std::clamp(desiredOutput[0], static_cast(Omin), static_cast(Omax)) : + desiredOutput[0]; + input = limitedOutput / dcGain; + fieldSet[0] = input - bias; + } + + if (stateCount > 0) { + if (std::abs(a.front()) < kMin_Res) { + if (std::abs(input) >= kMin_Res) { + throw InvalidParameterValue( + "transfer-function equilibrium requires a nonzero denominator constant"); + } + } else { + m_state[stateStart] = input * a.back() / a.front(); + } + for (index_t index = 1; index < stateCount; ++index) { + m_state[stateStart + index] = 0.0; + } + } + + m_state[rawOutputIndex] = rawOutput(input, m_state.data() + stateStart); + if (opFlags[USE_BLOCK_LIMITS]) { + GridBlock::rootCheck({input - bias}, + emptyStateData, + cLocalSolverMode, + CheckLevel::REVERSABLE_ONLY); + m_state[0] = vLimiter->clampOutput(m_state[rawOutputIndex]); + } + if (desiredOutput.empty()) { fieldSet[0] = m_state[0]; - prevInput = inputs[0] + bias; - } else { - m_state[0] = desiredOutput[0]; - // m_state[1] = (1.0 - m_T2 / m_T1) * desiredOutput[0] / K; - fieldSet[0] = desiredOutput[0] - bias; - prevInput = desiredOutput[0] / K; } + m_output = m_state[0]; + prevInput = input; } -// residual void TransferFunctionBlock::blockResidual(double input, - double didt, - const StateData& stateDataValue, - double resid[], + double /*didt*/, + const StateData& stateData, + double residual[], const SolverMode& sMode) { - auto loc = offsets.getLocations(stateDataValue, resid, sMode, this); - if (extraOutputState) { - } else { - for (size_t kk = 0; kk < a.size() - 1; ++kk) { - loc.destLoc[limiter_alg + kk] = - (-a[kk] * loc.diffStateLoc[kk]) + loc.diffStateLoc[kk + 1] + b[kk]; + const auto locations = offsets.getLocations(stateData, residual, sMode, this); + const index_t stateCount = order(); + if (hasDifferential(sMode)) { + for (index_t index = 0; index + 1 < stateCount; ++index) { + locations.destDiffLoc[index] = + locations.diffStateLoc[index + 1] - locations.dstateLoc[index]; + } + if (stateCount > 0) { + double finalDerivative = input + bias; + for (index_t index = 0; index < stateCount; ++index) { + finalDerivative -= (a[index] / a.back()) * locations.diffStateLoc[index]; + } + locations.destDiffLoc[stateCount - 1] = + finalDerivative - locations.dstateLoc[stateCount - 1]; + } + } + if (hasAlgebraic(sMode)) { + const index_t rawOutputIndex = limiter_alg; + locations.destLoc[rawOutputIndex] = + rawOutput(input + bias, locations.diffStateLoc) - locations.algStateLoc[rawOutputIndex]; + if (limiter_alg > 0) { + locations.destLoc[rawOutputIndex - 1] = + vLimiter->output(locations.algStateLoc[rawOutputIndex]) - + locations.algStateLoc[rawOutputIndex - 1]; } } - - // Loc.destLoc[limiter_alg] = Loc.diffStateLoc[limiter_diff] + m_T2 / m_T1 * - // (input + bias) - Loc.algStateLoc[limiter_alg]; - GridBlock::blockResidual(input, didt, stateDataValue, resid, sMode); } void TransferFunctionBlock::blockDerivative(double input, - double didt, - const StateData& stateDataValue, - double deriv[], + double /*didt*/, + const StateData& stateData, + double derivative[], const SolverMode& sMode) { - // auto offset = offsets.getDiffOffset (sMode); - // auto Aoffset = offsets.getAlgOffset (sMode); - // deriv[offset + limiter_diff] = K*(input + bias - sD.state[Aoffset + - // limiter_alg]) / m_T1; - if (opFlags[USE_RAMP_LIMITS]) { - GridBlock::blockDerivative(input, didt, stateDataValue, deriv, sMode); + if (!hasDifferential(sMode)) { + return; } + const auto locations = offsets.getLocations(stateData, derivative, sMode, this); + stateDerivative(input + bias, locations.diffStateLoc, locations.destDiffLoc); } -void TransferFunctionBlock::blockJacobianElements(double input, - double didt, - const StateData& stateDataValue, - MatrixData& matrixDataValue, - index_t argLoc, - const SolverMode& sMode) +void TransferFunctionBlock::blockAlgebraicUpdate(double input, + const StateData& stateData, + double update[], + const SolverMode& sMode) { - auto loc = offsets.getLocations(stateDataValue, sMode, this); - matrixDataValue.assign(loc.algOffset + 1, loc.algOffset + 1, -1); - - // md.assignCheck(Loc.algOffset + 1, argLoc, m_T2 / m_T1); - - GridBlock::blockJacobianElements(input, didt, stateDataValue, matrixDataValue, argLoc, sMode); - if (isAlgebraicOnly(sMode)) { + if (!hasAlgebraic(sMode)) { return; } - matrixDataValue.assign(loc.algOffset + 1, loc.diffOffset, 1); - // md.assign(arrayIndex, RowIndex, ColIndex, value) - - // md.assignCheck(Loc.diffOffset, argLoc, 1 / m_T1); - // md.assign(Loc.diffOffset, Loc.algOffset + 1, -1 / m_T1); - matrixDataValue.assign(loc.diffOffset, loc.diffOffset, -stateDataValue.cj); + const auto locations = offsets.getLocations(stateData, update, sMode, this); + const index_t rawOutputIndex = limiter_alg; + locations.destLoc[rawOutputIndex] = rawOutput(input + bias, locations.diffStateLoc); + if (limiter_alg > 0) { + locations.destLoc[rawOutputIndex - 1] = + vLimiter->output(locations.algStateLoc[rawOutputIndex]); + } } -double TransferFunctionBlock::step(CoreTime time, double inputA) +void TransferFunctionBlock::blockJacobianElements(double /*input*/, + double /*didt*/, + const StateData& stateData, + MatrixData& matrixData, + index_t inputLocation, + const SolverMode& sMode) { - const double timeDelta = time - prevTime; - double out; - const double input = inputA + bias; - // double ival, ival2; - if (timeDelta >= fabs(5.0)) { - m_state[2] = input; - } else if (timeDelta <= fabs(0.05)) { - // m_state[2] = m_state[2] + 1.0 / m_T1 * ((input + prevInput) / 2.0 - - // m_state[1]) * dt; - } else { - const double timeStep = 0.05; - double currentTime = prevTime + timeStep; - double intermediateInput = prevInput; - // double pin = prevInput; - // ival = m_state[2]; - // ival2 = m_state[1]; - while (currentTime < time) { - intermediateInput = intermediateInput + (((input - prevInput) / timeDelta) * timeStep); - // ival = ival + 1.0 / m_T1 * ((pin + in) / 2.0 - ival2) * tstep; - // ival2 = ival + m_T2 / m_T1 * (input); - currentTime += timeStep; - // pin = in; + const auto locations = offsets.getLocations(stateData, sMode, this); + const index_t stateCount = order(); + const double denominatorScale = a.back(); + if (hasDifferential(sMode)) { + for (index_t index = 0; index < stateCount; ++index) { + const index_t row = locations.diffOffset + index; + if (index + 1 < stateCount) { + matrixData.assign(row, row, -stateData.cj); + matrixData.assign(row, locations.diffOffset + index + 1, 1.0); + } else { + for (index_t column = 0; column < stateCount; ++column) { + matrixData.assign(row, + locations.diffOffset + column, + -a[column] / denominatorScale - + ((column == index) ? stateData.cj : 0.0)); + } + matrixData.assignCheckCol(row, inputLocation, 1.0); + } + } + } + if (hasAlgebraic(sMode)) { + const index_t rawOutputIndex = limiter_alg; + const index_t rawOutputLocation = locations.algOffset + rawOutputIndex; + matrixData.assign(rawOutputLocation, rawOutputLocation, -1.0); + if (stateCount == 0) { + matrixData.assignCheckCol(rawOutputLocation, + inputLocation, + K * b.front() / denominatorScale); + } else { + const double directTerm = b.back() / denominatorScale; + if (hasDifferential(sMode)) { + for (index_t index = 0; index < stateCount; ++index) { + matrixData.assign(rawOutputLocation, + locations.diffOffset + index, + K * + ((b[index] / denominatorScale) - + (directTerm * a[index] / denominatorScale))); + } + } + matrixData.assignCheckCol(rawOutputLocation, inputLocation, K * directTerm); + } + if (limiter_alg > 0) { + const index_t limitedOutputLocation = rawOutputLocation - 1; + matrixData.assign(limitedOutputLocation, limitedOutputLocation, -1.0); + matrixData.assign(limitedOutputLocation, rawOutputLocation, vLimiter->DoutDin()); } - // m_state[2] = ival + 1.0 / m_T1 * ((pin + input) / 2.0 - ival2) * (time - - // ct + tstep); } - // m_state[1] = m_state[2] + m_T2 / m_T1 * (input); +} - prevInput = input; - if (opFlags[HAS_LIMITS]) { - out = GridBlock::step(time, input); - } else { - out = K * m_state[1]; - m_state[0] = out; - prevTime = time; - m_output = out; +double TransferFunctionBlock::step(CoreTime time, double inputValue) +{ + const index_t stateCount = order(); + const index_t rawOutputIndex = limiter_alg; + const index_t stateStart = limiter_alg + 1; + const double input = inputValue + bias; + const double timeStep = time - prevTime; + if ((stateCount > 0) && (timeStep > 0.0)) { + const double halfStep = timeStep / 2.0; + std::vector systemMatrix(stateCount * stateCount, 0.0); + std::vector rightHandSide(stateCount, 0.0); + const double denominatorScale = a.back(); + for (index_t row = 0; row < stateCount; ++row) { + rightHandSide[row] = m_state[stateStart + row]; + for (index_t column = 0; column < stateCount; ++column) { + double systemValue = (row == column) ? 1.0 : 0.0; + double systemEntry = 0.0; + if (row + 1 < stateCount) { + systemEntry = (column == row + 1) ? 1.0 : 0.0; + } else { + systemEntry = -a[column] / denominatorScale; + } + systemValue -= halfStep * systemEntry; + rightHandSide[row] += halfStep * systemEntry * m_state[stateStart + column]; + systemMatrix[row * stateCount + column] = systemValue; + } + if (row + 1 == stateCount) { + rightHandSide[row] += halfStep * (prevInput + input); + } + } + for (index_t pivot = 0; pivot < stateCount; ++pivot) { + index_t pivotRow = pivot; + for (index_t row = pivot + 1; row < stateCount; ++row) { + if (std::abs(systemMatrix[row * stateCount + pivot]) > + std::abs(systemMatrix[pivotRow * stateCount + pivot])) { + pivotRow = row; + } + } + if (std::abs(systemMatrix[pivotRow * stateCount + pivot]) < kMin_Res) { + throw InvalidParameterValue("singular transfer-function timestep matrix"); + } + if (pivotRow != pivot) { + for (index_t column = pivot; column < stateCount; ++column) { + std::swap(systemMatrix[pivot * stateCount + column], + systemMatrix[pivotRow * stateCount + column]); + } + std::swap(rightHandSide[pivot], rightHandSide[pivotRow]); + } + const double pivotValue = systemMatrix[pivot * stateCount + pivot]; + for (index_t row = pivot + 1; row < stateCount; ++row) { + const double scale = systemMatrix[row * stateCount + pivot] / pivotValue; + for (index_t column = pivot; column < stateCount; ++column) { + systemMatrix[row * stateCount + column] -= + scale * systemMatrix[pivot * stateCount + column]; + } + rightHandSide[row] -= scale * rightHandSide[pivot]; + } + } + for (index_t row = stateCount; row-- > 0;) { + double value = rightHandSide[row]; + for (index_t column = row + 1; column < stateCount; ++column) { + value -= systemMatrix[row * stateCount + column] * m_state[stateStart + column]; + } + m_state[stateStart + row] = value / systemMatrix[row * stateCount + row]; + } + } + m_state[rawOutputIndex] = rawOutput(input, m_state.data() + stateStart); + if (opFlags[USE_BLOCK_LIMITS]) { + GridBlock::rootCheck({inputValue}, + emptyStateData, + cLocalSolverMode, + CheckLevel::REVERSABLE_ONLY); + m_state[0] = vLimiter->output(m_state[rawOutputIndex]); } - return out; + prevInput = input; + prevTime = time; + m_output = m_state[0]; + return m_output; } index_t TransferFunctionBlock::findIndex(std::string_view field, const SolverMode& sMode) const { - index_t ret = kInvalidLocation; - if (field == "m1") { - ret = offsets.getDiffOffset(sMode); - } else { - ret = GridBlock::findIndex(field, sMode); + if (field == "output") { + return offsets.getAlgOffset(sMode); } - return ret; + std::string prefix; + const int index = gmlc::utilities::stringOps::trailingStringInt(field, prefix, -1); + if ((index >= 0) && (prefix == "x")) { + if ((index >= 0) && (static_cast(index) < order())) { + return offsets.getDiffOffset(sMode) + static_cast(index); + } + } + return GridBlock::findIndex(field, sMode); } -// set parameters -void TransferFunctionBlock::set(std::string_view param, std::string_view val) +void TransferFunctionBlock::set(std::string_view param, std::string_view value) { if (param == "a") { - a = gmlc::utilities::str2vector(std::string{val}, 0); - } else if (param == "b") { - b = gmlc::utilities::str2vector(std::string{val}, 0); - } else { - GridBlock::set(param, val); + a = gmlc::utilities::str2vector(std::string{value}, 0.0); + if (a.empty()) { + throw InvalidParameterValue("transfer-function denominator"); + } + if (b.size() > a.size()) { + throw InvalidParameterValue("transfer-function numerator order"); + } + b.resize(a.size(), 0.0); + return; + } + if (param == "b") { + const auto numerator = gmlc::utilities::str2vector(std::string{value}, 0.0); + if (numerator.size() > a.size()) { + throw InvalidParameterValue("transfer-function numerator order"); + } + b = numerator; + b.resize(a.size(), 0.0); + return; } + GridBlock::set(param, value); } -void TransferFunctionBlock::set(std::string_view param, double val, units::unit unitType) +void TransferFunctionBlock::set(std::string_view param, double value, units::unit unitType) { - // param = GridDynSimulation::toLower(param); - std::string pstr; - const int num = gmlc::utilities::stringOps::trailingStringInt(param, pstr, -1); - if (pstr.length() == 1) { - switch (pstr[0]) { - case '#': - break; - case 'a': - case 't': - if (num >= 0) { - const auto index = static_cast(num); - if (index >= a.size()) { - a.resize(index + 1, 0); - b.resize(index + 1, 0); - } - a[index] = val; - } else { - throw(UnrecognizedParameter(param)); - } - break; - case 'b': - if (num >= 0) { - const auto index = static_cast(num); - if (index >= a.size()) { - a.resize(index + 1, 0); - b.resize(index + 1, 0); - } - b[index] = val; - } else { - throw(UnrecognizedParameter(param)); - } - break; - case 'k': - K = val; - break; - default: - throw(UnrecognizedParameter(param)); + std::string prefix; + const int coefficientIndex = gmlc::utilities::stringOps::trailingStringInt(param, prefix, -1); + if ((coefficientIndex >= 0) && ((prefix == "a") || (prefix == "b"))) { + const auto index = static_cast(coefficientIndex); + if (prefix == "a") { + if (index >= a.size()) { + a.resize(index + 1, 0.0); + b.resize(index + 1, 0.0); + } + a[index] = value; + } else { + if (index >= a.size()) { + throw InvalidParameterValue("transfer-function numerator order"); + } + b[index] = value; } + return; } - - if (param.empty() || param[0] == '#') { - // m_T1 = val; - } else { - GridBlock::set(param, val, unitType); - } + GridBlock::set(param, value, unitType); } -static stringVec gStateNames{"output", "Intermediate1", "intermediate2"}; - stringVec TransferFunctionBlock::localStateNames() const { - return gStateNames; + auto names = GridBlock::localStateNames(); + if (names.empty()) { + names.emplace_back("output"); + } else { + names.emplace_back("transfer_output"); + } + for (index_t index = 0; index < order(); ++index) { + names.emplace_back("x" + std::to_string(index)); + } + return names; } } // namespace griddyn::blocks diff --git a/src/griddyn/blocks/TransferFunctionBlock.h b/src/griddyn/blocks/TransferFunctionBlock.h index 3559d4e49..368704665 100644 --- a/src/griddyn/blocks/TransferFunctionBlock.h +++ b/src/griddyn/blocks/TransferFunctionBlock.h @@ -11,22 +11,35 @@ #include namespace griddyn::blocks { -/** @brief class implementing a generic transfer unction -block implementing \f$H(S)=\frac{K(b_0+b_1 s +/hdots +b_n s^n}{a_0+a_1 s +/hdots +a_n s^n}\f$ -it then converts it to observable canonical form as state space matrices for implementation as part -the solver - -*/ +/** + * @brief Proper single-input, single-output continuous-time transfer-function block. + * + * The coefficients are in ascending powers of @f$s@f$: + * @f[ + * G(s) = K\frac{b_0 + b_1s + \ldots + b_ns^n} + * {a_0 + a_1s + \ldots + a_ns^n}. + * @f] + * The numerator is zero-padded to the denominator order, so the block accepts a + * strictly proper transfer function as well as one with direct feedthrough. An + * improper transfer function is rejected at initialization. + * + * For a denominator of order @f$n@f$, GridDyn uses the controllable companion + * realization @f$\dot{x}_i=x_{i+1}@f$ for @f$i a; //!< lower time constant - std::vector b; //!< upper time constant - private: - // double rescale = 1; //!< containing the original $a_n$ for rescaling if - // coefficients are changed later - bool extraOutputState = false; //!< flag indicating that there is an extra state - //!< computation at the end due to direct dependence of B; + std::vector a; //!< denominator coefficients, ascending powers of s + std::vector b; //!< numerator coefficients, ascending powers of s public: /** constructor to add in the order of the transfer function @param[in] order the order of the transfer function @@ -63,6 +76,10 @@ are 0 const StateData& stateDataValue, double deriv[], const SolverMode& sMode) override; + virtual void blockAlgebraicUpdate(double input, + const StateData& stateDataValue, + double update[], + const SolverMode& sMode) override; virtual void blockResidual(double input, double didt, const StateData& stateDataValue, @@ -76,7 +93,22 @@ are 0 index_t argLoc, const SolverMode& sMode) override; virtual double step(CoreTime time, double inputA) override; - // virtual void setTime(CoreTime time){prevTime=time;}; virtual stringVec localStateNames() const override; + + /** + * @brief Evaluate the unbounded transfer function for state owned by a parent controller. + * + * This permits a composite controller to reuse the same companion realization + * while retaining ownership of its aggregate solver state and nonlinear limits. + */ + [[nodiscard]] double externalStateOutput(double input, const double state[]) const; + /** @brief Evaluate companion-state derivatives for state owned by a parent controller. */ + void externalStateDerivative(double input, const double state[], double derivative[]) const; + + private: + [[nodiscard]] index_t order() const; + void validateCoefficients(); + [[nodiscard]] double rawOutput(double input, const double state[]) const; + void stateDerivative(double input, const double state[], double derivative[]) const; }; } // namespace griddyn::blocks diff --git a/src/griddyn/blocks/ValueLimiter.h b/src/griddyn/blocks/ValueLimiter.h index 7e6c90427..6d578a3ff 100644 --- a/src/griddyn/blocks/ValueLimiter.h +++ b/src/griddyn/blocks/ValueLimiter.h @@ -9,12 +9,19 @@ #include namespace griddyn::blocks { -/** class that clamps a value between an upper and lower limit and maintains state of whether it - * is clamping or not +/** + * @brief Stateful hard value limiter used by GridBlock output limits. + * + * For requested value @f$v@f$, the static clamp is + * @f$\hat v=\min(\max(v,v_{min}),v_{max})@f$. @ref clampOutput evaluates + * that relation without state. The stateful @ref output, @ref deriv, and + * @ref DoutDin methods retain the active side of the limit so a DAE Jacobian + * has zero gain and zero derivative while clamped. `resetLevel` provides the + * release hysteresis used by root handling. */ class ValueLimiter { private: - double minVal = std::numeric_limits::min(); //!< Minimum value + double minVal = std::numeric_limits::lowest(); //!< Minimum value double maxVal = std::numeric_limits::max(); //!< maximum value double resetLevel = 0; //!< the amount the value has to go above or below a min or max to //!< be considered reset diff --git a/src/griddyn/blocks/blockLibrary.h b/src/griddyn/blocks/blockLibrary.h index 3df03d610..3c9b35c2d 100644 --- a/src/griddyn/blocks/blockLibrary.h +++ b/src/griddyn/blocks/blockLibrary.h @@ -18,6 +18,7 @@ #include "FunctionBlock.h" #include "IntegralBlock.h" #include "LutBlock.h" +#include "LeadLagBlock.h" #include "NullBlock.h" #include "PidBlock.h" #include "TransferFunctionBlock.h" diff --git a/src/griddyn/governors/GovernorTgov1.cpp b/src/griddyn/governors/GovernorTgov1.cpp index f5bde0c0b..e4db3f778 100644 --- a/src/griddyn/governors/GovernorTgov1.cpp +++ b/src/griddyn/governors/GovernorTgov1.cpp @@ -13,7 +13,6 @@ #include "utilities/MatrixData.hpp" #include #include -#include #include namespace griddyn::governors { @@ -57,6 +56,10 @@ void GovernorTgov1::dynObjectInitializeA(CoreTime time0, std::uint32_t flags) (T1 <= 0.0) || (T3 <= 0.0) || (Pmax < Pmin)) { throw InvalidParameterValue("TGOV1 R, T1, T3, or valve limits"); } + // ANDES exposes the turbine output as its state and uses + // T3 * pmdot = valve - pm - T2 * valvedot. This is the output-state + // realization of (1 - T2 s)/(1 + T3 s), hence the negative lead time. + turbineTransfer.setParameters(-T2, T3); GovernorIeeeSimple::dynObjectInitializeA(time0, flags); } @@ -108,8 +111,9 @@ void GovernorTgov1::derivative(const IOdata& inputs, (-governorState[1] + inputs[govpSetInLocation] - (K * (omega - 1.0))) / T1; } - loc.destDiffLoc[0] = - (loc.diffStateLoc[1] - loc.diffStateLoc[0] - (T2 * loc.destDiffLoc[1])) / T3; + loc.destDiffLoc[0] = turbineTransfer.outputStateDerivative(governorState[1], + governorState[0], + loc.destDiffLoc[1]); } void GovernorTgov1::timestep(CoreTime time, const IOdata& inputs, const SolverMode& /*sMode*/) @@ -159,8 +163,12 @@ if (opFlags.test (uses_deadband)) if (opFlags[POWER_LIMITED]) { matrixData.assign(referenceIndex + 1, referenceIndex + 1, -stateData.cj); - matrixData.assign(referenceIndex, referenceIndex, -(1 / T3) - stateData.cj); - matrixData.assign(referenceIndex, referenceIndex + 1, 1 / T3); + matrixData.assign(referenceIndex, + referenceIndex, + turbineTransfer.derivativeStateJacobian() - stateData.cj); + matrixData.assign(referenceIndex, + referenceIndex + 1, + turbineTransfer.derivativeInputJacobian()); } else { matrixData.assignCheckCol(referenceIndex + 1, inputLocs[govpSetInLocation], 1 / T1); matrixData.assign(referenceIndex + 1, referenceIndex + 1, -(1 / T1) - stateData.cj); @@ -168,12 +176,21 @@ if (opFlags.test (uses_deadband)) matrixData.assign(referenceIndex + 1, inputLocs[govOmegaInLocation], -K / T1); } - matrixData.assign(referenceIndex, referenceIndex + 1, (1 + (T2 / T1)) / T3); - matrixData.assignCheckCol(referenceIndex, inputLocs[govpSetInLocation], -T2 / T1 / T3); + matrixData.assign(referenceIndex, + referenceIndex + 1, + turbineTransfer.derivativeInputJacobian() - + (turbineTransfer.outputStateInputDerivativeJacobian() / T1)); + matrixData.assignCheckCol(referenceIndex, + inputLocs[govpSetInLocation], + turbineTransfer.outputStateInputDerivativeJacobian() / T1); if (linkOmega) { - matrixData.assign(referenceIndex, inputLocs[govOmegaInLocation], K * T2 / T1 / T3); + matrixData.assign(referenceIndex, + inputLocs[govOmegaInLocation], + -K * turbineTransfer.outputStateInputDerivativeJacobian() / T1); } - matrixData.assign(referenceIndex, referenceIndex, -(1 / T3) - stateData.cj); + matrixData.assign(referenceIndex, + referenceIndex, + turbineTransfer.derivativeStateJacobian() - stateData.cj); } } diff --git a/src/griddyn/governors/GovernorTgov1.h b/src/griddyn/governors/GovernorTgov1.h index 6573c36e5..5f8ed43f7 100644 --- a/src/griddyn/governors/GovernorTgov1.h +++ b/src/griddyn/governors/GovernorTgov1.h @@ -7,6 +7,7 @@ #pragma once #include "GovernorIeeeSimple.h" +#include "blocks/LeadLag.h" #include #include @@ -15,9 +16,14 @@ namespace griddyn::governors { * @brief PSS/e TGOV1 turbine governor. * * The valve state follows \f$\dot v=(-v+P_{ref}-\Delta\omega/R)/T_1\f$, - * with the PSS/e VMIN/VMAX limits. The turbine output is the lead-lag - * \f$(1+sT_2)/(1+sT_3)\f$ response of the valve state, minus - * \f$D_t\Delta\omega\f$. + * with the PSS/e VMIN/VMAX limits. ANDES publishes the turbine section in + * its output-state form, + * \f[ + * T_3\dot p_m=v-p_m-T_2\dot v, + * \f] + * which has transfer function \f$(1-sT_2)/(1+sT_3)\f$ from valve position + * \f$v\f$ to turbine output \f$p_m\f$. The mechanical output is + * \f$p_m-D_t\Delta\omega\f$. * * These equations match ANDES v2.0.0 * `andes/models/governor/tgov1.py`, TGOV1Model. @@ -26,6 +32,9 @@ class GovernorTgov1: public GovernorIeeeSimple { public: protected: double Dt = 0.0; //!< speed damping constant + /** Turbine lead-lag kernel; this governor retains ownership of its output-state realization. */ + blocks::LeadLagKernel turbineTransfer; + public: explicit GovernorTgov1(const std::string& objName = "govTgov1_#"); virtual CoreObject* clone(CoreObject* obj = nullptr) const override; diff --git a/test/componentTests/testBlocks.cpp b/test/componentTests/testBlocks.cpp index b6f2983fa..56a49b5c3 100644 --- a/test/componentTests/testBlocks.cpp +++ b/test/componentTests/testBlocks.cpp @@ -5,10 +5,13 @@ */ #include "../gtestHelper.h" +#include "core/CoreExceptions.h" #include "core/ObjectFactory.hpp" #include "gmlc/utilities/TimeSeriesMulti.hpp" #include "gmlc/utilities/vectorOps.hpp" #include "griddyn/Relay.h" +#include "griddyn/blocks/RampLimiter.h" +#include "griddyn/blocks/ValueLimiter.h" #include "griddyn/blocks/blockLibrary.h" #include "griddyn/simulation/Diagnostics.h" #include @@ -209,6 +212,124 @@ TEST_F(BlockTests, DeadbandBlockTest) EXPECT_EQ(ret, 0); } +TEST_F(BlockTests, TransferFunctionInitializationAndStep) +{ + // G(s) = 1 / (1 + s). The companion state equals the output because + // the numerator has no direct feedthrough. + TransferFunctionBlock block({1.0, 1.0}, {1.0}); + block.dynInitializeA(0.0, 0U); + IOdata fieldSet; + block.dynInitializeB({2.0}, {}, fieldSet); + EXPECT_NEAR(fieldSet[0], 2.0, 1e-12); + EXPECT_NEAR(block.getBlockOutput(), 2.0, 1e-12); + + // The local stepping path uses an implicit trapezoidal update and treats + // successive inputs as a linear interpolation. The input therefore + // ramps from 2 to 0 over this first 0.01 s interval. + const double output = block.step(0.01, 0.0); + EXPECT_NEAR(output, 2.0 / 1.005, 1e-12); + + // Desired-output initialization back-solves through the finite DC gain. + TransferFunctionBlock initializedBlock({2.0, 1.0}, {4.0}); + initializedBlock.dynInitializeA(0.0, 0U); + initializedBlock.dynInitializeB({}, {3.0}, fieldSet); + EXPECT_NEAR(fieldSet[0], 1.5, 1e-12); + EXPECT_NEAR(initializedBlock.getBlockOutput(), 3.0, 1e-12); +} + +TEST_F(BlockTests, TransferFunctionOutputLimitAndParameterValidation) +{ + TransferFunctionBlock limitedBlock({1.0, 1.0}, {1.0}); + limitedBlock.set("max", 0.5); + limitedBlock.set("min", -0.5); + limitedBlock.dynInitializeA(0.0, 0U); + IOdata fieldSet; + limitedBlock.dynInitializeB({2.0}, {}, fieldSet); + EXPECT_NEAR(limitedBlock.getBlockOutput(), 0.5, 1e-12); + // The output limiter releases when the unbounded companion output returns + // within range; it does not freeze the transfer-function state. + EXPECT_NEAR(limitedBlock.step(10.0, 0.0), 1.0 / 3.0, 1e-12); + + TransferFunctionBlock invalidBlock; + EXPECT_THROW(invalidBlock.set("b", "1,2,3"), InvalidParameterValue); + invalidBlock.set("a", "1,0"); + EXPECT_THROW(invalidBlock.dynInitializeA(0.0, 0U), InvalidParameterValue); +} + +TEST_F(BlockTests, LeadLagKernelAndBlock) +{ + // G(s) = 2(1 + 0.5s)/(1 + 2s). The equilibrium state equals the input. + LeadLagKernel kernel(0.5, 2.0, 2.0); + EXPECT_TRUE(kernel.isValid()); + EXPECT_NEAR(kernel.output(3.0, 1.0), 3.0, 1e-12); + EXPECT_NEAR(kernel.derivative(3.0, 1.0), 1.0, 1e-12); + EXPECT_NEAR(kernel.outputInputJacobian(), 0.5, 1e-12); + EXPECT_NEAR(kernel.outputStateJacobian(), 1.5, 1e-12); + + LeadLagBlock block(2.0, 0.5, 2.0); + block.dynInitializeA(0.0, 0U); + IOdata fieldSet; + block.dynInitializeB({3.0}, {}, fieldSet); + EXPECT_NEAR(fieldSet[0], 6.0, 1e-12); + EXPECT_NEAR(block.getBlockOutput(), 6.0, 1e-12); + + // The lag state has the exact zero-order-hold response in the local-step path. + EXPECT_NEAR(block.step(2.0, 1.0), 2.0 + (3.0 / std::exp(1.0)), 1e-12); + + LeadLagBlock invalidBlock; + EXPECT_THROW(invalidBlock.set("tb", 0.0), InvalidParameterValue); +} + +TEST_F(BlockTests, LutBlockInterpolationInitializationAndLimits) +{ + LutBlock block; + block.set("lut", "0,0;1,2;2,2"); + EXPECT_NEAR(block.computeValue(-1.0), 0.0, 1e-12); + EXPECT_NEAR(block.computeValue(0.5), 1.0, 1e-12); + EXPECT_NEAR(block.computeValue(1.5), 2.0, 1e-12); + EXPECT_NEAR(block.computeValue(3.0), 2.0, 1e-12); + + block.dynInitializeA(0.0, 0U); + IOdata fieldSet; + block.dynInitializeB({0.5}, {}, fieldSet); + EXPECT_NEAR(fieldSet[0], 1.0, 1e-12); + EXPECT_NEAR(block.getBlockOutput(), 1.0, 1e-12); + + LutBlock initializedBlock; + initializedBlock.set("lut", "0,0;1,2;2,4"); + initializedBlock.dynInitializeA(0.0, 0U); + initializedBlock.dynInitializeB({}, {3.0}, fieldSet); + EXPECT_NEAR(fieldSet[0], 1.5, 1e-12); + EXPECT_NEAR(initializedBlock.getBlockOutput(), 3.0, 1e-12); + + LutBlock invalidBlock; + EXPECT_THROW(invalidBlock.set("lut", "0,0,1"), InvalidParameterValue); + EXPECT_THROW(invalidBlock.set("lut", "0,0;0,1"), InvalidParameterValue); +} + +TEST_F(BlockTests, FunctionBlockRejectsUnknownFunction) +{ + EXPECT_THROW(FunctionBlock("not_a_function"), InvalidParameterValue); + + FunctionBlock block("sin"); + EXPECT_THROW(block.set("func", "not_a_function"), InvalidParameterValue); + block.dynInitializeA(0.0, 0U); + IOdata fieldSet; + block.dynInitializeB({0.5}, {}, fieldSet); + EXPECT_NEAR(block.getBlockOutput(), std::sin(0.5), 1e-12); +} + +TEST_F(BlockTests, DefaultLimitersDoNotClampNegativeValues) +{ + ValueLimiter valueLimiter; + EXPECT_NEAR(valueLimiter.clampOutput(-1.0), -1.0, 1e-12); + valueLimiter.changeLimitActivation(-1.0); + EXPECT_FALSE(valueLimiter.isActive()); + + RampLimiter rampLimiter; + EXPECT_NEAR(rampLimiter.clampOutputRamp(-1.0), -1.0, 1e-12); +} + using blockdescpair = std::pair>>; std::vector makeBlockParameterMap() @@ -227,8 +348,11 @@ std::vector makeBlockParameterMap() std::make_pair("d", 0.28), std::make_pair("t", 0.2)}}, {"control", {std::make_pair("t1", 0.2), std::make_pair("t2", 0.1)}}, + {"leadlag", {std::make_pair("tb", 0.2), std::make_pair("ta", 0.1)}}, {"function", {std::make_pair("gain", kPI), std::make_pair("bias", -0.05)}}, {"func", {std::make_pair("arg", 2.35)}}, + {"tf", {}}, + {"lut", {}}, }; } @@ -239,6 +363,8 @@ std::map>> {"function", {std::make_pair("func", "sin")}}, {"func", {std::make_pair("func", "pow")}}, {"db", {std::make_pair("flags", "shifted")}}, + {"tf", {{"a", "1,3,2"}, {"b", "1"}}}, + {"lut", {{"lut", "-1,-2;1,2;2,3"}}}, }; } @@ -321,7 +447,7 @@ TEST_P(BlockCompareTests, CompareBlockTest) EXPECT_EQ(ret, 0); } -INSTANTIATE_TEST_SUITE_P(AllBlocks, BlockCompareTests, ::testing::Range(0, 11)); +INSTANTIATE_TEST_SUITE_P(AllBlocks, BlockCompareTests, ::testing::Range(0, 14)); #ifdef GRIDDYN_ENABLE_CVODE /** test the control block if they can handle a differential only Jacobian and an algebraic only From 9f459070541064388cd3574f2b052daf865b2234 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sat, 29 Aug 2026 00:01:34 +0000 Subject: [PATCH 2/6] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/griddyn/blocks/blockLibrary.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/griddyn/blocks/blockLibrary.h b/src/griddyn/blocks/blockLibrary.h index 3c9b35c2d..81a621db1 100644 --- a/src/griddyn/blocks/blockLibrary.h +++ b/src/griddyn/blocks/blockLibrary.h @@ -17,8 +17,8 @@ #include "FilteredDerivativeBlock.h" #include "FunctionBlock.h" #include "IntegralBlock.h" -#include "LutBlock.h" #include "LeadLagBlock.h" +#include "LutBlock.h" #include "NullBlock.h" #include "PidBlock.h" #include "TransferFunctionBlock.h" From 111dd658097120ddd9d13484aa06efc3005427fa Mon Sep 17 00:00:00 2001 From: Philip Top Date: Fri, 28 Aug 2026 19:12:35 -0700 Subject: [PATCH 3/6] Potential fix for pull request finding 'CodeQL / Multiplication result converted to larger type' Co-authored-by: Copilot Autofix powered by AI <62310815+github-advanced-security[bot]@users.noreply.github.com> --- src/griddyn/blocks/TransferFunctionBlock.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/griddyn/blocks/TransferFunctionBlock.cpp b/src/griddyn/blocks/TransferFunctionBlock.cpp index 00363ce9b..4af2e0b21 100644 --- a/src/griddyn/blocks/TransferFunctionBlock.cpp +++ b/src/griddyn/blocks/TransferFunctionBlock.cpp @@ -334,7 +334,7 @@ double TransferFunctionBlock::step(CoreTime time, double inputValue) const double timeStep = time - prevTime; if ((stateCount > 0) && (timeStep > 0.0)) { const double halfStep = timeStep / 2.0; - std::vector systemMatrix(stateCount * stateCount, 0.0); + std::vector systemMatrix(static_cast(stateCount) * static_cast(stateCount), 0.0); std::vector rightHandSide(stateCount, 0.0); const double denominatorScale = a.back(); for (index_t row = 0; row < stateCount; ++row) { From 0907e8745c5440730845aef1d969c6f17933d1c7 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sat, 29 Aug 2026 02:13:15 +0000 Subject: [PATCH 4/6] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/griddyn/blocks/TransferFunctionBlock.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/griddyn/blocks/TransferFunctionBlock.cpp b/src/griddyn/blocks/TransferFunctionBlock.cpp index 4af2e0b21..18aeec21e 100644 --- a/src/griddyn/blocks/TransferFunctionBlock.cpp +++ b/src/griddyn/blocks/TransferFunctionBlock.cpp @@ -334,7 +334,9 @@ double TransferFunctionBlock::step(CoreTime time, double inputValue) const double timeStep = time - prevTime; if ((stateCount > 0) && (timeStep > 0.0)) { const double halfStep = timeStep / 2.0; - std::vector systemMatrix(static_cast(stateCount) * static_cast(stateCount), 0.0); + std::vector systemMatrix(static_cast(stateCount) * + static_cast(stateCount), + 0.0); std::vector rightHandSide(stateCount, 0.0); const double denominatorScale = a.back(); for (index_t row = 0; row < stateCount; ++row) { From dd951e97ba422a57bc1a7cf58f2b9281524f62d0 Mon Sep 17 00:00:00 2001 From: Philip Top Date: Fri, 28 Aug 2026 19:22:42 -0700 Subject: [PATCH 5/6] fix warnings --- src/griddyn/blocks/LeadLag.h | 12 ++++--- src/griddyn/blocks/LutBlock.cpp | 6 ++-- src/griddyn/blocks/TransferFunctionBlock.cpp | 34 ++++++++++++-------- src/griddyn/blocks/TransferFunctionBlock.h | 4 +-- src/griddyn/governors/GovernorTgov1.cpp | 1 + 5 files changed, 34 insertions(+), 23 deletions(-) diff --git a/src/griddyn/blocks/LeadLag.h b/src/griddyn/blocks/LeadLag.h index 55249f618..506560cb7 100644 --- a/src/griddyn/blocks/LeadLag.h +++ b/src/griddyn/blocks/LeadLag.h @@ -105,13 +105,17 @@ class LeadLagKernel { { return (input - outputState + (Ta * inputDerivative)) / Tb; } - /** @return @f$\partial y/\partial u@f$. */ + /** @brief Return the output partial derivative with respect to input, @f$\partial y/\partial + * u@f$. */ [[nodiscard]] double outputInputJacobian() const { return K * Ta / Tb; } - /** @return @f$\partial y/\partial x@f$. */ + /** @brief Return the output partial derivative with respect to lag state, @f$\partial + * y/\partial x@f$. */ [[nodiscard]] double outputStateJacobian() const { return K * (1.0 - (Ta / Tb)); } - /** @return @f$\partial\dot{x}/\partial u@f$. */ + /** @brief Return the lag-state derivative partial with respect to input, + * @f$\partial\dot{x}/\partial u@f$. */ [[nodiscard]] double derivativeInputJacobian() const { return 1.0 / Tb; } - /** @return @f$\partial\dot{x}/\partial x@f$. */ + /** @brief Return the lag-state derivative partial with respect to state, + * @f$\partial\dot{x}/\partial x@f$. */ [[nodiscard]] double derivativeStateJacobian() const { return -1.0 / Tb; } /** @return the coefficient of @f$\dot u@f$ in the output-state realization. */ [[nodiscard]] double outputStateInputDerivativeJacobian() const { return Ta / Tb; } diff --git a/src/griddyn/blocks/LutBlock.cpp b/src/griddyn/blocks/LutBlock.cpp index 27835a5d5..ca4c03e90 100644 --- a/src/griddyn/blocks/LutBlock.cpp +++ b/src/griddyn/blocks/LutBlock.cpp @@ -198,10 +198,10 @@ LutBlock::LookupResult LutBlock::evaluate(double input) const { validateTable(lut); if (lut.size() == 1 || input <= lut.front().first) { - return {lut.front().second, 0.0}; + return {.value = lut.front().second, .slope = 0.0}; } if (input >= lut.back().first) { - return {lut.back().second, 0.0}; + return {.value = lut.back().second, .slope = 0.0}; } const auto upper = @@ -210,7 +210,7 @@ LutBlock::LookupResult LutBlock::evaluate(double input) const }); const auto lower = std::prev(upper); const double slope = (upper->second - lower->second) / (upper->first - lower->first); - return {lower->second + ((input - lower->first) * slope), slope}; + return {.value = lower->second + ((input - lower->first) * slope), .slope = slope}; } double LutBlock::inverseValue(double value) const diff --git a/src/griddyn/blocks/TransferFunctionBlock.cpp b/src/griddyn/blocks/TransferFunctionBlock.cpp index 00363ce9b..a1e8c52f3 100644 --- a/src/griddyn/blocks/TransferFunctionBlock.cpp +++ b/src/griddyn/blocks/TransferFunctionBlock.cpp @@ -289,7 +289,7 @@ void TransferFunctionBlock::blockJacobianElements(double /*input*/, for (index_t column = 0; column < stateCount; ++column) { matrixData.assign(row, locations.diffOffset + column, - -a[column] / denominatorScale - + (-a[column] / denominatorScale) - ((column == index) ? stateData.cj : 0.0)); } matrixData.assignCheckCol(row, inputLocation, 1.0); @@ -334,7 +334,13 @@ double TransferFunctionBlock::step(CoreTime time, double inputValue) const double timeStep = time - prevTime; if ((stateCount > 0) && (timeStep > 0.0)) { const double halfStep = timeStep / 2.0; - std::vector systemMatrix(stateCount * stateCount, 0.0); + const auto matrixIndex = [stateCount](index_t row, index_t column) { + return (static_cast(row) * static_cast(stateCount)) + + static_cast(column); + }; + std::vector systemMatrix(static_cast(stateCount) * + static_cast(stateCount), + 0.0); std::vector rightHandSide(stateCount, 0.0); const double denominatorScale = a.back(); for (index_t row = 0; row < stateCount; ++row) { @@ -349,7 +355,7 @@ double TransferFunctionBlock::step(CoreTime time, double inputValue) } systemValue -= halfStep * systemEntry; rightHandSide[row] += halfStep * systemEntry * m_state[stateStart + column]; - systemMatrix[row * stateCount + column] = systemValue; + systemMatrix[matrixIndex(row, column)] = systemValue; } if (row + 1 == stateCount) { rightHandSide[row] += halfStep * (prevInput + input); @@ -358,27 +364,27 @@ double TransferFunctionBlock::step(CoreTime time, double inputValue) for (index_t pivot = 0; pivot < stateCount; ++pivot) { index_t pivotRow = pivot; for (index_t row = pivot + 1; row < stateCount; ++row) { - if (std::abs(systemMatrix[row * stateCount + pivot]) > - std::abs(systemMatrix[pivotRow * stateCount + pivot])) { + if (std::abs(systemMatrix[matrixIndex(row, pivot)]) > + std::abs(systemMatrix[matrixIndex(pivotRow, pivot)])) { pivotRow = row; } } - if (std::abs(systemMatrix[pivotRow * stateCount + pivot]) < kMin_Res) { + if (std::abs(systemMatrix[matrixIndex(pivotRow, pivot)]) < kMin_Res) { throw InvalidParameterValue("singular transfer-function timestep matrix"); } if (pivotRow != pivot) { for (index_t column = pivot; column < stateCount; ++column) { - std::swap(systemMatrix[pivot * stateCount + column], - systemMatrix[pivotRow * stateCount + column]); + std::swap(systemMatrix[matrixIndex(pivot, column)], + systemMatrix[matrixIndex(pivotRow, column)]); } std::swap(rightHandSide[pivot], rightHandSide[pivotRow]); } - const double pivotValue = systemMatrix[pivot * stateCount + pivot]; + const double pivotValue = systemMatrix[matrixIndex(pivot, pivot)]; for (index_t row = pivot + 1; row < stateCount; ++row) { - const double scale = systemMatrix[row * stateCount + pivot] / pivotValue; + const double scale = systemMatrix[matrixIndex(row, pivot)] / pivotValue; for (index_t column = pivot; column < stateCount; ++column) { - systemMatrix[row * stateCount + column] -= - scale * systemMatrix[pivot * stateCount + column]; + systemMatrix[matrixIndex(row, column)] -= + scale * systemMatrix[matrixIndex(pivot, column)]; } rightHandSide[row] -= scale * rightHandSide[pivot]; } @@ -386,9 +392,9 @@ double TransferFunctionBlock::step(CoreTime time, double inputValue) for (index_t row = stateCount; row-- > 0;) { double value = rightHandSide[row]; for (index_t column = row + 1; column < stateCount; ++column) { - value -= systemMatrix[row * stateCount + column] * m_state[stateStart + column]; + value -= systemMatrix[matrixIndex(row, column)] * m_state[stateStart + column]; } - m_state[stateStart + row] = value / systemMatrix[row * stateCount + row]; + m_state[stateStart + row] = value / systemMatrix[matrixIndex(row, row)]; } } m_state[rawOutputIndex] = rawOutput(input, m_state.data() + stateStart); diff --git a/src/griddyn/blocks/TransferFunctionBlock.h b/src/griddyn/blocks/TransferFunctionBlock.h index 368704665..fdf0c6f17 100644 --- a/src/griddyn/blocks/TransferFunctionBlock.h +++ b/src/griddyn/blocks/TransferFunctionBlock.h @@ -90,9 +90,9 @@ are 0 double didt, const StateData& stateDataValue, MatrixData& matrixDataValue, - index_t argLoc, + index_t inputLocation, const SolverMode& sMode) override; - virtual double step(CoreTime time, double inputA) override; + virtual double step(CoreTime time, double inputValue) override; virtual stringVec localStateNames() const override; /** diff --git a/src/griddyn/governors/GovernorTgov1.cpp b/src/griddyn/governors/GovernorTgov1.cpp index e4db3f778..6ac2a9936 100644 --- a/src/griddyn/governors/GovernorTgov1.cpp +++ b/src/griddyn/governors/GovernorTgov1.cpp @@ -13,6 +13,7 @@ #include "utilities/MatrixData.hpp" #include #include +#include #include namespace griddyn::governors { From 43dcd115b5b06fe714ef989a93fa35c7ee9147b5 Mon Sep 17 00:00:00 2001 From: Philip Top Date: Sat, 29 Aug 2026 04:32:14 -0700 Subject: [PATCH 6/6] cpplint fixes --- src/griddyn/blocks/LeadLagBlock.cpp | 1 + src/griddyn/blocks/LutBlock.cpp | 1 + 2 files changed, 2 insertions(+) diff --git a/src/griddyn/blocks/LeadLagBlock.cpp b/src/griddyn/blocks/LeadLagBlock.cpp index bfad0ce09..105cafdfb 100644 --- a/src/griddyn/blocks/LeadLagBlock.cpp +++ b/src/griddyn/blocks/LeadLagBlock.cpp @@ -11,6 +11,7 @@ #include "core/CoreObjectTemplates.hpp" #include "utilities/MatrixData.hpp" #include +#include namespace griddyn::blocks { LeadLagBlock::LeadLagBlock(const std::string& objName): GridBlock(objName) diff --git a/src/griddyn/blocks/LutBlock.cpp b/src/griddyn/blocks/LutBlock.cpp index ca4c03e90..a95f8ca1e 100644 --- a/src/griddyn/blocks/LutBlock.cpp +++ b/src/griddyn/blocks/LutBlock.cpp @@ -16,6 +16,7 @@ #include #include #include +#include namespace griddyn::blocks { LutBlock::LutBlock(const std::string& objName): GridBlock(objName)