API — ICoreBlocks/ICoreMath/ControlSystems
The public contract of 19 header(s) under src/ICoreBlocks/ICoreMath/ControlSystems — 19 class/struct definition(s), 121 declaration(s). Each section shows the header's banner and its public (and protected-virtual) surface exactly as the file writes it.
ICoreIIREmulator.h#
src/ICoreBlocks/ICoreMath/ControlSystems/ICoreIIREmulator.h
ICoreIIREmulator#
ICoreIIREmulator.h:9 · class · pImpl · 12 declaration(s)
class ICoreIIREmulator {
public:
// explicit ICoreIIREmulator(const std::vector<double>& numCoefficients, const std::vector<double>& denCoefficients);
explicit ICoreIIREmulator(const ICoreTransferFunction& tf);
double step(const double& uk);
void reset();
void setInitialConditions(const std::vector<double>& new_u_ic, const std::vector<double>& new_y_ic);
// Seed the filter from a DIRECT-FORM II delay line, newest first -- the vector
// Simulink's Discrete Transfer Fcn and Discrete Filter call "Initial states"
// (both report FilterStructure = "Direct form II"). This emulator is direct
// form I: it holds past INPUTS and past OUTPUTS, which are different
// quantities, so the two cannot simply be assigned to each other.
//
// The conversion, with w[-m] the m-th newest given state and 0 outside the
// vector, and the coefficients already normalized by a0:
//
// u[-m] = w[-m] + SUM_j a_j * w[-m-j]
// y[-m] = SUM_j b_j * w[-m-j]
//
// MEASURED against R2026a rather than derived from the documentation: for
// num [0.5 -0.2 0.3], den [1 -0.7 0.25] and InitialStates [1 0], Simulink's
// zero-input response begins 0.15, 0.28, and so does this seeding. The two
// realizations agree sample for sample, which is the only claim worth making
// -- an initial STATE is realization-specific and could not be copied across.
//
// Also measured: "InitialDenominatorStates" moves nothing on either block, so
// there is no second vector to carry.
void seedFromDirectFormIIStates(const std::vector<double>& states);
// The same conversion as data, for the code generators: they bake the seeded
// histories into the exported filter and never call step(). Index i is the
// (i+1)-th newest sample, matching what step() expects to find.
[[nodiscard]] std::vector<double> seededInputHistory() const;
[[nodiscard]] std::vector<double> seededOutputHistory() const;
// Emulators are stored BY VALUE in std::vector inside three block headers
// (Discrete Transfer Function, Transfer Function, Zero-Pole), so the type
// has to stay both copyable and nothrow-movable: a vector only grows by
// moving when the move cannot throw, and the implicit copy is deleted by
// the unique_ptr residue. All four are written out in the .cpp.
ICoreIIREmulator(const ICoreIIREmulator& other);
ICoreIIREmulator& operator=(const ICoreIIREmulator& other);
ICoreIIREmulator(ICoreIIREmulator&& other) noexcept;
ICoreIIREmulator& operator=(ICoreIIREmulator&& other) noexcept;
~ICoreIIREmulator();
private:
class Impl; // the two-line residue; state lives here
std::unique_ptr<Impl> impl;
};ICoreStateSpace.h#
src/ICoreBlocks/ICoreMath/ControlSystems/ICoreStateSpace.h
///////////////////////////// / Constructors ////////////////////////////
ICoreStateSpace#
ICoreStateSpace.h:11 · class · pImpl · 38 declaration(s)
class ICoreStateSpace {
public:
///////////////////////////////
/// Constructors
//////////////////////////////
ICoreStateSpace();
explicit ICoreStateSpace(const ICoreMatrix& A_candidate, const ICoreMatrix& B_candidate, const ICoreMatrix& C_candidate, const ICoreMatrix& D_candidate, const double& Ts_candidate = -1.0);
explicit ICoreStateSpace(const ICoreMatrix& A_candidate, const ICoreMatrix& Bu_candidate, const ICoreMatrix& Bf_candidate,
const ICoreMatrix& C_candidate, const ICoreMatrix& Du_candidate, const ICoreMatrix& Df_candidate, const double& Ts_candidate = -1.0);
static ICoreStateSpace fromTransferFunction(const ICoreTransferFunction& tf);
///////////////////////////////
/// Matrices Assignment
//////////////////////////////
bool tryUpdatingMatrices(const ICoreMatrix &A_candidate, const ICoreMatrix &B_candidate,
const ICoreMatrix &C_candidate, const ICoreMatrix &D_candidate,
const double &Ts_candidate);
bool tryUpdatingMatrices(const ICoreMatrix &A_candidate, const ICoreMatrix &Bu_candidate,
const ICoreMatrix &Bf_candidate, const ICoreMatrix &C_candidate,
const ICoreMatrix &Du_candidate, const ICoreMatrix &Df_candidate,
const double &Ts_candidate);
///////////////////////////////
/// Dynamics
//////////////////////////////
[[nodiscard]] ICoreMatrix computeStateEvolution_noMatrixSizeCheck(const ICoreMatrix& x, const ICoreMatrix& u) const;
[[nodiscard]] ICoreMatrix computeStateEvolution_noMatrixSizeCheck(const ICoreMatrix &x, const ICoreMatrix &u, const ICoreMatrix &f) const;
[[nodiscard]] ICoreMatrix computeOutput_noMatrixSizeCheck(const ICoreMatrix& x, const ICoreMatrix& u) const;
[[nodiscard]] ICoreMatrix computeOutput_noMatrixSizeCheck(const ICoreMatrix& x, const ICoreMatrix& u, const ICoreMatrix &f) const;
///////////////////////////////
/// Analysis
//////////////////////////////
[[nodiscard]] bool checkDimensions(const ICoreMatrix& x, const ICoreMatrix& u) const;
bool checkDimensions_nonlinear(const ICoreMatrix &x, const ICoreMatrix &u, const ICoreMatrix &f) const;
[[nodiscard]] std::vector<ICoreComplexVariable> getEigenValues() const;
[[nodiscard]] ICoreMatrix getControllabilityMatrix() const;
[[nodiscard]] bool isControllable() const;
[[nodiscard]] ICoreMatrix getObservabilityMatrix() const;
[[nodiscard]] bool isObservable() const;
///////////////////////////////
/// Getters
//////////////////////////////
[[nodiscard]] ICoreMatrix getA() const;
[[nodiscard]] ICoreMatrix getB() const;
[[nodiscard]] ICoreMatrix getBu() const;
[[nodiscard]] ICoreMatrix getBf() const;
[[nodiscard]] ICoreMatrix getCombinedBuBf() const;
[[nodiscard]] ICoreMatrix getC() const;
[[nodiscard]] ICoreMatrix getD() const;
[[nodiscard]] ICoreMatrix getDu() const;
[[nodiscard]] ICoreMatrix getDf() const;
[[nodiscard]] ICoreMatrix getCombinedDuDf() const;
[[nodiscard]] size_t getOrder() const;
[[nodiscard]] size_t getNumberOfInputs() const;
[[nodiscard]] size_t get_m_u() const;
[[nodiscard]] size_t get_m_f() const;
[[nodiscard]] size_t getNumberOfOutputs() const;
[[nodiscard]] double getSamplingTime() const;
[[nodiscard]] double getTs() const;
[[nodiscard]] bool isContinuous() const;
void print() const;
// A state space is passed and stored by value all over the model layer, so
// it stays copyable. The unique_ptr residue below deletes the implicit copy
// operations, so all four are written out in the .cpp.
ICoreStateSpace(const ICoreStateSpace& other);
ICoreStateSpace& operator=(const ICoreStateSpace& other);
ICoreStateSpace(ICoreStateSpace&& other) noexcept;
ICoreStateSpace& operator=(ICoreStateSpace&& other) noexcept;
~ICoreStateSpace();
private:
class Impl; // the two-line residue; state lives here
std::unique_ptr<Impl> impl;
};ICoreTransferFunction.h#
src/ICoreBlocks/ICoreMath/ControlSystems/ICoreTransferFunction.h
MATLAB's rule for the sample time two models share when they meet (toolbox board T0.2 (1), T1.1): a STATIC GAIN -- numerator and denominator both constants -- adopts the other model's sample time, and otherwise the two must agree, or the combination is refused with MATLAB's own sentence ("Sampling times must agree."). The three operators above and feedback() take their Ts from here, so a scalar on the LEFT of a discrete model no longer hands the product a continuous sample time. Answers false with
whyNotset when the two disagree; the operators then fall back to the left operand's Ts, which is why a caller that can refuse asks first.
ICoreTransferFunction#
ICoreTransferFunction.h:14 · class · pImpl · 34 declaration(s)
class ICoreTransferFunction {
public:
ICoreTransferFunction();
ICoreTransferFunction(const ICorePolynomial &num_candidate, const ICorePolynomial &den_candidate, const double& Ts_candidate = -1);
static ICoreTransferFunction fromPolesZeros(const std::vector<ICoreComplexVariable>& zeros,
const std::vector<ICoreComplexVariable>& poles,
const double& gain = 1.0, const double& Ts = -1);
static ICoreTransferFunction fromSISOStateSpace(const ICoreStateSpace &stateSpace);
static ICoreTransferFunction fromSISOStateSpace(const ICoreMatrix& A, const ICoreMatrix& B,
const ICoreMatrix& C, const ICoreMatrix& D, double Ts = -1);
bool tryUpdatingCoefficients(const ICorePolynomial &num_candidate, const ICorePolynomial &den_candidate,
const double &Ts_candidate = -1);
ICoreTransferFunction operator*(const ICoreTransferFunction &other) const;
ICoreTransferFunction operator+(const ICoreTransferFunction &other) const;
ICoreTransferFunction operator-(const ICoreTransferFunction &other) const;
// MATLAB's rule for the sample time two models share when they meet
// (toolbox board T0.2 (1), T1.1): a STATIC GAIN -- numerator
// and denominator both constants -- adopts the other model's sample time,
// and otherwise the two must agree, or the combination is refused with
// MATLAB's own sentence ("Sampling times must agree."). The three
// operators above and feedback() take their Ts from here, so a scalar
// on the LEFT of a discrete model no longer hands the product a
// continuous sample time. Answers false with `whyNot` set when the two
// disagree; the operators then fall back to the left operand's Ts, which
// is why a caller that can refuse asks first.
static bool sampleTimesAgree(const ICoreTransferFunction& a, const ICoreTransferFunction& b,
double& Ts, std::string* whyNot = nullptr);
// Numerator degree at most the denominator's. An IMPROPER model is a
// value since T1.1 -- `tf('s')` is one, and so is every polynomial in s
// built from it -- and the operations MATLAB refuses for one (a time
// response, a zoh/foh/impulse discretization, an explicit state-space
// realization) ask this before they start.
[[nodiscard]] bool isProper() const;
// Numerator and denominator both constants -- MATLAB's "Static gain."
[[nodiscard]] bool isStaticGain() const;
// Every numerator coefficient zero: the model that displays as `0`.
[[nodiscard]] bool hasZeroNumerator() const;
// den/num with the same sample time, which is what `1/G` and `G^-1` mean
// in MATLAB. Nothing is cancelled and nothing is normalized, so
// inverse().inverse() is the model it started from coefficient for
// coefficient.
[[nodiscard]] ICoreTransferFunction inverse() const;
// G^n for an integer n: a repeated product for n > 0, the static gain 1
// (with this model's sample time) for n == 0, and inverse()^|n| below
// zero -- MATLAB's mpower on a SISO model, which multiplies rather than
// cancels, so (s+1)/(s+1) squared is a fourth-order fraction on both
// sides.
[[nodiscard]] ICoreTransferFunction power(int exponent) const;
// G/(1 + G*H). MATLAB's default and the classical negative loop.
[[nodiscard]] ICoreTransferFunction feedback(const ICoreTransferFunction &H) const;
// The same loop with the sign of the junction written out:
// `sign` is -1 for the negative loop above and +1 for MATLAB's
// `feedback(G, H, +1)`, which closes G/(1 - G*H). The two differ only
// in that one sign, so they are ONE formula with a parameter rather
// than two -- a second copy is where the positive loop eventually
// stops agreeing with the negative one about a shared numerator
// (toolbox board T1.6).
[[nodiscard]] ICoreTransferFunction feedback(const ICoreTransferFunction &H,
const double &sign) const;
void poleZeroCancellation(const double& tol);
void normalize();
[[nodiscard]] ICoreComplexVariable evaluate(const ICoreComplexVariable &s) const;
[[nodiscard]] std::vector<ICoreComplexVariable> zeros() const;
[[nodiscard]] std::vector<ICoreComplexVariable> poles() const;
[[nodiscard]] bool isStable() const;
[[nodiscard]] double getDcGain() const;
[[nodiscard]] ICoreComplexVariable frequencyResponse(const double& omega) const;
[[nodiscard]] std::vector<ICoreBodePoint> bode(const double& w_start, const double& w_end, const int& numOfPoints) const;
[[nodiscard]] std::vector<ICoreComplexVariable> nyquist(const double& w_start, const double& w_end, const int& numOfPoints) const;
[[nodiscard]] ICorePolynomial getNumerator() const;
[[nodiscard]] ICorePolynomial getDenominator() const;
[[nodiscard]] double getSamplingTime() const;
[[nodiscard]] double getTs() const;
// ---- the PID facet (toolbox board T1.30) -----------------
//
// ⚠ A `pid` IS A TRANSFER FUNCTION HERE, and that is the row's own
// recommendation rather than a shortcut: T1.30 asked whether the console
// should grow a `Kind::Pid` and answered no -- a `tf` carrying a DISPLAY
// flag and the four readable gains, converting to a plain `tf` on any
// arithmetic. MATLAB does the same thing at the value level (`pid * tf`
// is a `tf` there too), so the two agree on what a controller BECOMES the
// moment it is combined with anything.
//
// The facet is data, not behaviour: nothing in this class reads it, the
// three operators build fresh models that do not carry it, and the only
// consumers are the display (ICoreControlSystemsFormatting) and the
// console's postfix member reads (`C.Kp`).
//
// ⚠ IT IS ALSO PART OF THE STORAGE. A console variable is stored as its
// canonical text, so `toCanonical` writes `pid(...)`/`pidstd(...)` for a
// model carrying the facet and `parseTransferFunction` reads it back --
// otherwise a pid would silently become a `tf` the moment it was assigned
// to a name.
enum class PidForm {
None, // an ordinary transfer function
Parallel, // MATLAB's `pid`: Kp + Ki/s + Kd*s/(Tf*s + 1)
Standard // MATLAB's `pidstd`: Kp*(1 + 1/(Ti*s) + Td*s/((Td/N)*s + 1))
};
// The four gains, in MATLAB's own order for the form: Kp, Ki, Kd, Tf for
// the parallel one and Kp, Ti, Td, N for the standard one.
void setPidForm(const PidForm& form, const double& first, const double& second,
const double& third, const double& fourth);
[[nodiscard]] PidForm pidForm() const;
[[nodiscard]] double pidGain(const size_t& index) const; // 0..3
[[nodiscard]] int getOrder() const;
void print() const;
// Transfer functions are handed around and stored by value, so the type
// stays copyable; the unique_ptr residue below deletes the implicit copy
// operations, so all four are written out in the .cpp.
ICoreTransferFunction(const ICoreTransferFunction& other);
ICoreTransferFunction& operator=(const ICoreTransferFunction& other);
ICoreTransferFunction(ICoreTransferFunction&& other) noexcept;
ICoreTransferFunction& operator=(ICoreTransferFunction&& other) noexcept;
~ICoreTransferFunction();
// The Pade approximation of a pure delay e^{-tau s} (toolbox row T1.32),
// as MATLAB's `[num, den] = pade(tau, n)` answers it: DESCENDING
// coefficient rows, and monic -- MATLAB scales both by the leading
// denominator coefficient, so the pair is exact rational arithmetic and
// not a fit.
//
// den(s) = sum_k ((2n-k)! n!) / ((2n)! k! (n-k)!) (tau s)^k
// num(s) = den(-s)
//
// so the numerator is the denominator with alternating signs and nothing
// else -- which is why they are built once here and why `pade(tau, 0)` is
// the pair (1, 1): the delay is dropped, not approximated.
//
// Answers false with `whyNot` set for a negative delay or order, which is
// what MATLAB refuses too ("Value must be nonnegative").
static bool padeCoefficients(const double& delay, const double& order,
std::vector<double>& numeratorDescending,
std::vector<double>& denominatorDescending,
std::string* whyNot = nullptr);
private:
class Impl; // the two-line residue; state lives here
std::unique_ptr<Impl> impl;
};ICoreZpk.h#
src/ICoreBlocks/ICoreMath/ControlSystems/ICoreZpk.h
ICoreZpk#
ICoreZpk.h:30 · class · pImpl · nested Factor · 13 declaration(s)
MATLAB's zpk model as a console VALUE (toolbox board T1.3, T0.2 (5), X3): the zeros, the poles, the scalar gain and the sample time, held as themselves rather than as the coefficients they multip...
class ICoreZpk {
public:
ICoreZpk();
ICoreZpk(const std::vector<ICoreComplexVariable>& zeros,
const std::vector<ICoreComplexVariable>& poles,
const double& gain = 1.0, const double& Ts = -1.0);
// The roots of a transfer function, with the gain MATLAB's own zpk(sys)
// answers: the leading numerator coefficient over the leading denominator
// one. A model whose numerator is identically zero is the zpk with no
// zeros and gain 0, which is what MATLAB answers too.
static ICoreZpk fromTransferFunction(const ICoreTransferFunction& tf);
[[nodiscard]] ICoreTransferFunction toTransferFunction() const;
[[nodiscard]] std::vector<ICoreComplexVariable> getZeros() const;
[[nodiscard]] std::vector<ICoreComplexVariable> getPoles() const;
[[nodiscard]] double getGain() const;
[[nodiscard]] double getTs() const;
// No zeros and no poles: MATLAB's "Static gain." display, and the value
// that adopts the other operand's sample time when two models meet
// (T0.2 (1), ICoreTransferFunction::sampleTimesAgree).
[[nodiscard]] bool isStaticGain() const;
// Zeros and poles GROUPED for display: equal roots collapsed to a
// multiplicity, and each conjugate pair collapsed to the real quadratic
// MATLAB prints instead of two complex linear factors. `tolerance` is
// relative to the root's own magnitude, which is what makes -1 +/- 1e-9i
// print as (s+1)^2 the way MATLAB prints it.
struct Factor {
// A linear factor (s - (root + b i)) when `quadratic` is false, and
// the real quadratic (s^2 + b s + c) when it is true.
//
// ⚠ `b` MEANS TWO THINGS, one per arm, and the arm is `quadratic`.
// A linear factor's `b` is the IMAGINARY part of its root, which is
// zero for every factor a conjugate-paired model produces and
// non-zero only for the unpaired complex root MATLAB warns about and
// then prints as (s-(1+2i)). Folding it into `root` would need a
// complex member in a struct two display functions read, and giving
// the unpaired case its own struct would duplicate the multiplicity
// logic for a case that is one line of display.
bool quadratic = false;
double root = 0.0; // linear: the real part of the root
double b = 0.0, c = 0.0; // see above
size_t multiplicity = 1;
};
static std::vector<Factor> groupForDisplay(const std::vector<ICoreComplexVariable>& roots,
const double& tolerance = 1e-8);
ICoreZpk(const ICoreZpk& other);
ICoreZpk& operator=(const ICoreZpk& other);
ICoreZpk(ICoreZpk&& other) noexcept;
ICoreZpk& operator=(ICoreZpk&& other) noexcept;
~ICoreZpk();
private:
class Impl; // the two-line residue; state lives here
std::unique_ptr<Impl> impl;
};ICoreBodePoint.h#
src/ICoreBlocks/ICoreMath/ControlSystems/HelperObjects/ICoreBodePoint.h
bode() returns a std::vector of these, so the type has to stay copyable. A unique_ptr member deletes the implicit copy operations, so both are written out in the .cpp -- see the note there on what the residue costs.
ICoreBodePoint#
ICoreBodePoint.h:7 · class · pImpl · 9 declaration(s)
class ICoreBodePoint {
public:
ICoreBodePoint(const double& omega,
const double& magnitude,
const double& magnitude_dB,
const double& phase
);
// bode() returns a std::vector of these, so the type has to stay copyable.
// A unique_ptr member deletes the implicit copy operations, so both are
// written out in the .cpp -- see the note there on what the residue costs.
ICoreBodePoint(const ICoreBodePoint& other);
ICoreBodePoint& operator=(const ICoreBodePoint& other);
ICoreBodePoint(ICoreBodePoint&& other) noexcept;
ICoreBodePoint& operator=(ICoreBodePoint&& other) noexcept;
[[nodiscard]] double getOmega() const;
[[nodiscard]] double getMagnitude() const;
[[nodiscard]] double getMagnitudeDB() const;
[[nodiscard]] double getPhase() const;
~ICoreBodePoint();
private:
class Impl; // the two-line residue; state lives here
std::unique_ptr<Impl> impl;
};ICoreBlocksMerging.h#
src/ICoreBlocks/ICoreMath/ControlSystems/StaticMethods/ICoreBlocksMerging.h
MATLAB's
append(sys1, sys2)(toolbox row T1.34): the two systems side by side and not connected at all -- block-diagonal A, B and C, D. The result has every input and every output of both, in order, which is what makes it the building block the index forms ofseriesandparallelpick channels out of.⚠ It is NOT
mergeInParallel, which SUMS two systems driven by one input. Both put the A matrices on a diagonal and there the resemblance stops: parallel gives one input and one output, append gives two of each, and confusing them answers a model of the right ORDER with the wrong shape.
ICoreBlocksMerging#
ICoreBlocksMerging.h:6 · class · 5 declaration(s)
class ICoreBlocksMerging {
public:
static ICoreStateSpace mergeInSeries(const ICoreStateSpace& G1, const ICoreStateSpace& G2);
static ICoreStateSpace mergeInParallel(const ICoreStateSpace& G1, const ICoreStateSpace& G2);
// MATLAB's `append(sys1, sys2)` (toolbox row T1.34): the two systems side
// by side and not connected at all -- block-diagonal A, B and C, D. The
// result has every input and every output of both, in order, which is what
// makes it the building block the index forms of `series` and `parallel`
// pick channels out of.
//
// ⚠ It is NOT `mergeInParallel`, which SUMS two systems driven by one
// input. Both put the A matrices on a diagonal and there the resemblance
// stops: parallel gives one input and one output, append gives two of
// each, and confusing them answers a model of the right ORDER with the
// wrong shape.
static ICoreStateSpace appendDiagonal(const ICoreStateSpace& G1, const ICoreStateSpace& G2);
// Closed-loop realization for a purely additive summing junction: e = r + H(y), y = G(e). (For
// the classical negative-feedback junction e = r - H(y), bake the sign into H's own C/D — the
// formula here does not assume either sign.) Resolves the algebraic loop created by G and H's
// direct feedthrough terms; returns false (leaving outClosedLoop untouched) if that loop is
// unresolvable, i.e. 1 - D_G*D_H is (numerically) singular.
static bool mergeInFeedback(const ICoreStateSpace& G, const ICoreStateSpace& H, ICoreStateSpace& outClosedLoop);
// Closed-loop realization for an "empty"/unity feedback path: e = referenceSign*r +
// feedbackSign*y, y = G(e) (G's own output wired straight back into the summing junction, with
// no feedback block in between). referenceSign/feedbackSign are each +1 or -1, matching which
// port of the summing junction the reference and the loop-closing connection land on (see
// ICoreStateSpaceReductionHelpers::combinerInputSign) — e.g. classical negative unity feedback
// is referenceSign=+1, feedbackSign=-1. Resolves the algebraic loop created by G's own direct
// feedthrough term; returns false (leaving outClosedLoop untouched) if that loop is
// unresolvable, i.e. 1 - feedbackSign*D_G is (numerically) singular.
static bool mergeInUnityFeedback(const ICoreStateSpace& G, double referenceSign, double feedbackSign, ICoreStateSpace& outClosedLoop);
};
};ICoreControlSignals.h#
src/ICoreBlocks/ICoreMath/ControlSystems/StaticMethods/ICoreControlSignals.h
ICoreControlSignals#
ICoreControlSignals.h:18 · class · 1 declaration(s)
The Control System Toolbox test-input generator (the toolbox board's T1.35): the periodic signal lsim is driven with.
class ICoreControlSignals {
public:
// Which of MATLAB's three classes of signal to generate. MATLAB selects on
// the first TWO characters of the word, so "si", "sin", "sine" all name the
// sine and "sq"/"square" the square wave -- the console reads it the same
// way, which is why this enum is produced by gensigType() below rather than
// by an exact-match table.
enum class Kind { Sine, Square, Pulse };
// Reads MATLAB's TYPE word. Answers false, with `whyNot` set, for a word
// shorter than two characters or one whose first two match none of the
// three -- MATLAB's own precondition (`length(type)>1`) and its own
// switch, in that order.
static bool gensigType(const std::string& word, Kind& kindOut, std::string* whyNot);
// [u, t] = gensig(type, tau, Tf, Ts). `tau` is the period in seconds,
// `duration` the length of the signal (MATLAB's default 5*tau) and
// `sampleTime` the spacing of the samples (MATLAB's default tau/64). Both
// answers are COLUMNS of the same height, and the height is that of
// MATLAB's `(0:Ts:Tf)'` -- floor(Tf/Ts) + 1 samples, not a rounded count.
//
// The three signals, exactly MATLAB's:
// Sine sin(2*pi*t/tau)
// Square 1 where rem(t, tau) >= tau/2, else 0
// Pulse 1 where rem(t, tau) < (1 - 1000*eps)*Ts, else 0
//
// Answers false with `whyNot` set when tau or Ts is not positive or the
// duration is negative -- where MATLAB divides by zero and answers a
// signal of NaNs, or builds an empty column, without saying anything.
static bool gensig(const Kind& kind, const double& tau, const double& duration,
const double& sampleTime, ICoreMatrix& signalOut, ICoreMatrix& timeOut,
std::string* whyNot);
};
};ICoreMatlabFrequencyResponse.h#
src/ICoreBlocks/ICoreMath/ControlSystems/StaticMethods/ICoreMatlabFrequencyResponse.h
ICoreMatlabFrequencyResponse#
ICoreMatlabFrequencyResponse.h:58 · class · nested Spec, Result · 0 declaration(s)
MATLAB's bode, nyquist, nichols and sigma on MATLAB's own frequency grid (toolbox board rows T1.18 and T1.19).
class ICoreMatlabFrequencyResponse {
public:
// MATLAB's GRADE, spelled by the verb that passes it (`freqresp.m`:
// "GRADE should be set to 1 for NYQUIST, 2 for NICHOLS, and 3 for BODE";
// `sigmaresp.m` passes 4).
enum class Verb { Nyquist = 1, Nichols = 2, Bode = 3, Sigma = 4 };
// What the caller asked for, which is MATLAB's `wspec`: nothing at all,
// a RANGE `{wmin, wmax}` the grid is built inside, or a GRID of the
// caller's own. Only the last escapes the transcription above -- and it is
// the only one MATLAB answers in the order it was given.
struct Spec {
enum class Kind { Automatic, Range, Grid };
Kind kind = Kind::Automatic;
double minimum = 0.0; // Range: wmin, 0 for "unspecified"
double maximum = 0.0; // Range: wmax, infinite for "none"
std::vector<double> frequencies; // Grid
};
// Everything the four verbs read, at every frequency of the grid. Callers
// take the two or three columns their verb answers: `bode` magnitude and
// phase, `nyquist` real and imaginary, `nichols` the same pair as `bode`,
// `sigma` magnitude alone (the one singular value of a SISO model).
//
// `phase` is in RADIANS and is CONTINUOUS -- MATLAB's `bode` converts it
// to degrees at the very end (`rad2deg`), and nothing unwraps it after the
// fact: it is accumulated per pole and per zero, so a phase that has run
// past -180 degrees stays there.
struct Result {
std::vector<double> frequencies;
std::vector<double> magnitude; // absolute, NOT dB
std::vector<double> phase; // radians, continuous
std::vector<double> real;
std::vector<double> imaginary;
};
static bool compute(const ICoreTransferFunction& model, const Verb& verb, const Spec& spec,
const bool& focusFromMagnitude, Result& out, std::string& whyNot);
};
};ICoreMatlabResponse.h#
src/ICoreBlocks/ICoreMath/ControlSystems/StaticMethods/ICoreMatlabResponse.h
ICoreMatlabResponse#
ICoreMatlabResponse.h:44 · class · nested Grid, Characteristics · 1 declaration(s)
MATLAB's step, impulse and initial on MATLAB's own time grid (toolbox board rows T1.13, T1.14, T1.16).
class ICoreMatlabResponse {
public:
enum class Input { Step, Impulse, Initial };
// `dt` and `Tf` for one response, and how many samples the pair implies.
struct Grid {
double dt = 0.0;
double finalTime = 0.0;
size_t samples = 0; // including t = 0
bool settlingIsMatlabs = true; // false when MATLAB would extend it
};
// timegrid.m. `initialState` is the x0 of an `initial` response and is
// ignored by the other two. `finalTime` <= 0 asks for the automatic
// horizon, where the answer's `settlingIsMatlabs` is false.
//
// ⚠ A DISCRETE model runs timegrid TOO, and MATLAB's does not. MATLAB
// steps a discrete model at its own Ts and leaves the horizon entirely to
// the compiled settling detector (trespSetUp's `isCT && isempty(dt)`
// guard), which left this console answering NO SAMPLES for
// `step(c2d(sys, Ts))` -- a refusal, not a response, and the reason T1.17
// could not read a discrete step at all. The horizon now comes from
// timegrid's own rule read in the z plane: the continuous pole of a
// discrete one is s = log(z)/Ts, and every line of timegrid then applies
// unchanged over those s. Only the COUNT is taken from it; the SPACING is
// always the model's own Ts. See the note in the .cpp.
static Grid grid(const ICoreStateSpace& sys, const Input& input,
const ICoreMatrix& initialState, const double& finalTime);
// The grid a caller's own time VECTOR implies. MATLAB rebuilds a uniform
// grid from `t(1)`, `t(2) - t(1)` and `t(end)` rather than using the
// vector as given (getTimeInfo.m + timeresp.m), so a non-uniform vector is
// refused on both sides and a uniform one may still come back a sample
// longer or shorter than it went in.
static bool gridFromTimeVector(const ICoreMatrix& times, Grid& out, std::string& whyNot);
// The response itself, over `grid`. `y` is `samples x 1` for a SISO model
// and `x` the state trajectory `samples x nx`; both are answered whether
// or not the caller wants x, because they cost one simulation.
//
// The discretization is MATLAB's, and it is NOT one map for the three:
// step and initial are ZOH (`utDiscretizeZOH`), impulse is IMPULSE
// INVARIANT (`utDiscretizeIMP`) -- "a discrete model whose response to a
// unit-area pulse of length Ts and height 1/Ts matches the continuous-time
// response to the Dirac impulse", which for a delay-free model is
// y[k] = C e^{A k Ts} B with the state read POST-impulse.
// MATLAB's `lsim(sys, u, t[, x0])` (T1.15). `u` and `times` are Ns-long
// vectors; `x0` may be empty. `y` comes back Ns x 1 and `x` Ns x nx, the
// CONTINUOUS state trajectory.
//
// ⚠ MATLAB PICKS THE INTERPOLATION RULE FROM THE INPUT, and the choice is
// a HARD SWITCH rather than a blend: `utSelectInterp.m`, fourteen readable
// lines, normalises by `range = max(u) - min(u)` and calls the input
// SMOOTH -- so FOH -- when `max(abs(diff(u)))` is at most 0.75 of that
// range, and jumpy (ZOH) otherwise. Measured on tf(1,[1 2 1]) at
// Ts = 0.1: sin(t) matches a FOH simulation to 0.0 and a ZOH one to
// 1.7e-2; a unit step at t = 1 matches ZOH to 0.0 and FOH to 1.8e-2. So
// it picks one per input and gets it exactly.
static bool linearSimulation(const ICoreStateSpace& sys, const ICoreMatrix& u,
const ICoreMatrix& times, const ICoreMatrix& initialState,
ICoreMatrix& y, ICoreMatrix& x, std::string& whyNot);
// The same simulation with the interpolation rule NAMED instead of picked:
// `intersample` is "zoh" or "foh", and nothing else is accepted.
//
// ⚠ It exists because `utSelectInterp` is a rule for a USER's `lsim`, and
// the System Identification Toolbox never lets it run: every filter inside
// `inival_time.m` passes its rule explicitly, and the two are not the same
// choice. `svf` filters the OUTPUT through 'foh' and the INPUT through the
// record's own intersample; `srivc` filters the output through 'foh' and
// the input through 'zoh' whatever the record says. A binary input picked
// ZOH and a smooth one FOH would answer a different initial model for the
// same record (T6.2), so the estimator asks for the rule by name.
static bool linearSimulation(const ICoreStateSpace& sys, const ICoreMatrix& u,
const ICoreMatrix& times, const ICoreMatrix& initialState,
const std::string& intersample,
ICoreMatrix& y, ICoreMatrix& x, std::string& whyNot);
static bool simulate(const ICoreStateSpace& sys, const Input& input,
const ICoreMatrix& initialState, const Grid& grid,
ICoreMatrix& times, ICoreMatrix& y, ICoreMatrix& x,
std::string& whyNot);
// MATLAB's `stepinfo` (toolbox row T1.17): the nine characteristics read
// off a step response, in the ORDER MATLAB answers them.
//
// It lives beside the grid because that is where it lives in MATLAB --
// `stepinfo.m` sits in the same `controllib/engine` directory as
// `timegrid.m` -- and because the two halves of the row are the same
// question asked twice: `stepinfo(y, t)` reads a response someone else
// computed, `stepinfo(sys)` reads the one grid() and simulate() compute
// right here. Sharing the file is what keeps them from drifting.
//
// ⚠ THE FIELD ORDER IS NOT THE ORDER THE DOCUMENTATION LISTS, and a
// record is an ORDERED list, so getting it wrong prints a struct MATLAB
// never prints. `fieldnames(stepinfo(sys))` on R2026a is
// RiseTime, TransientTime, SettlingTime, SettlingMin, SettlingMax,
// Overshoot, Undershoot, Peak, PeakTime -- TransientTime is SECOND.
//
// ⚠ `undershoot` IS NEGATIVE ZERO for every response that never dips
// below its initial value, because MATLAB computes it as
// `-100 * min(0, min(yrel))` and the C++ product carries the sign bit the
// same way. It is not a defect to be tidied away: a record that answers
// +0 there disagrees with MATLAB in the sign bit, which `1/x` can see.
struct Characteristics {
double riseTime = 0.0;
double transientTime = 0.0;
double settlingTime = 0.0;
double settlingMin = 0.0;
double settlingMax = 0.0;
double overshoot = 0.0;
double undershoot = 0.0;
double peak = 0.0;
double peakTime = 0.0;
// stepinfo.m's undocumented SECOND output. Computed because the rise
// time is their difference and SettlingMin/SettlingMax are read from
// t >= riseTimeHigh; not answered to the console (see the catalog
// note on `stepinfo`).
double riseTimeLow = 0.0;
double riseTimeHigh = 0.0;
};
// LocalGetInfo, transcribed line for line from
// `controllib/engine/stepinfo.m` (R2026a, P. Gahinet, 1986-2021).
//
// `sampleTime` is stepinfo's undocumented 'Ts' option and it does ONE
// thing: a nonzero value turns off every linear interpolation, so a
// discrete response reports the sample it crossed at rather than a time
// between two samples. `@DynamicSystem/stepinfo.m` passes `abs(sys.Ts)`
// for exactly that reason, which is why the model form takes it too.
//
// ⚠ TWO TOLERANCES WITH ONE OPTION NAME, and they are different numbers.
// TransientTime settles against `threshold * max|y - yfinal|` and
// SettlingTime against `threshold * |yfinal - yinit|`; the two coincide
// only when the largest error is the initial one, which is exactly the
// case a monotone response satisfies and an overshooting one does not.
static bool characteristics(const ICoreMatrix& y, const ICoreMatrix& times,
const double& yFinal, const double& yInit,
const double& settlingThreshold, const double& riseLow,
const double& riseHigh, const double& sampleTime,
Characteristics& out, std::string& whyNot);
};
};ICoreModelArithmetic.h#
src/ICoreBlocks/ICoreMath/ControlSystems/StaticMethods/ICoreModelArithmetic.h
ICoreModelArithmetic#
ICoreModelArithmetic.h:28 · class · 3 declaration(s)
MATLAB's ARITHMETIC on LTI models -- G1 * G2, G1 + G2, G1 - G2, G1 / G2, inv(G), G' and G^n (toolbox board row T1.5).
class ICoreModelArithmetic {
public:
// sys1 * sys2 -- series with the LEFT factor upstream of the output:
// A = [A1, B1*C2; 0, A2] B = [B1*D2; B2]
// C = [C1, D1*C2] D = D1*D2
// `whyNot` is set (and the answer left default) when the inner dimensions
// disagree or the sample times do.
static bool multiply(const ICoreStateSpace& left, const ICoreStateSpace& right,
ICoreStateSpace& out, std::string& whyNot);
// sys1 + sys2 and sys1 - sys2 -- one input driving both, the outputs
// summed: A = blkdiag(A1, A2), B = [B1; B2], C = [C1, +/-C2],
// D = D1 +/- D2.
static bool addOrSubtract(const ICoreStateSpace& a, const ICoreStateSpace& b,
const bool& subtract, ICoreStateSpace& out, std::string& whyNot);
// inv(sys). MATLAB answers a DESCRIPTOR model here -- `inv(ss(-3,1,2,0.5))`
// is a two-state model with E = diag(1, 0), measured on R2026a -- and this
// console's ICoreStateSpace has no E matrix. So the answer is the ordinary
// explicit inverse, which is the SAME MODEL in a smaller realization:
// A - B*inv(D)*C, B*inv(D), -inv(D)*C, inv(D)
// It exists only when D is square and invertible; a model with singular
// feedthrough (every strictly proper one) is refused by name rather than
// answered as something else.
static bool inverse(const ICoreStateSpace& sys, ICoreStateSpace& out, std::string& whyNot);
// sys1 / sys2, i.e. sys1 * inv(sys2), with inverse()'s restriction.
static bool divide(const ICoreStateSpace& a, const ICoreStateSpace& b,
ICoreStateSpace& out, std::string& whyNot);
// sys^n for an integer n: repeated multiply above zero, the static
// identity gain at zero, and inverse() first below it.
static bool power(const ICoreStateSpace& sys, const int& exponent,
ICoreStateSpace& out, std::string& whyNot);
// sys' -- the conjugate transpose. For a CONTINUOUS model it is
// (-A', -C', B', D'), measured on R2026a, which is G(-s)^T. A DISCRETE
// one is G(1/z)^T, and MATLAB realizes it as a descriptor with a singular
// E (measured: a 2-state discrete model transposes to a 3-state E-model),
// so it is refused here for inverse()'s reason.
static bool conjugateTranspose(const ICoreStateSpace& sys,
ICoreStateSpace& out, std::string& whyNot);
// sys.' -- the PLAIN transpose, which is (A', C', B', D') and needs no
// descriptor on either time base, because it transposes the I/O map
// without substituting the variable. For a SISO model it is the model
// itself, and a SISO transfer function's `.'` is likewise itself.
static ICoreStateSpace plainTranspose(const ICoreStateSpace& sys);
// G' for a SISO transfer function: G(-s) when continuous -- the sign of
// every odd power flipped -- and G(1/z) when discrete, which is both
// coefficient rows padded to one length and reversed. `G.'` is G itself
// for a SISO model and does not need a function.
static ICoreTransferFunction conjugateTranspose(const ICoreTransferFunction& model);
};
};ICoreModelFrequencyReadings.h#
src/ICoreBlocks/ICoreMath/ControlSystems/StaticMethods/ICoreModelFrequencyReadings.h
ICoreModelFrequencyReadings#
ICoreModelFrequencyReadings.h:40 · class · 0 declaration(s)
The frequency-domain readings of a model, toolbox row T1.21 -- freqresp, evalfr, bandwidth, getPeakGain and the H-infinity half of norm.
class ICoreModelFrequencyReadings {
public:
// G at each w. Continuous: G(jw). Discrete: G(e^{jwTs}).
static bool frequencyResponse(const ICoreTransferFunction& model,
const std::vector<double>& frequencies,
std::vector<double>& real, std::vector<double>& imaginary,
std::string* whyNot = nullptr);
// G at a RAW point of the complex plane -- no jw and no e^{jwTs}.
static void evaluateAt(const ICoreTransferFunction& model, const double& pointReal,
const double& pointImaginary, double& real, double& imaginary);
// The peak of |G| along the frequency axis, and the frequency it is at.
// Continuous: exact, through the stationary points of |G(jw)|^2 as a
// rational function of w^2. Discrete: the same question over
// theta in [0, pi], refined to machine precision.
static bool peakGain(const ICoreTransferFunction& model, double& peak, double& frequency,
std::string* whyNot = nullptr);
// The first frequency at which |G| has fallen `dropDb` below its value at
// zero frequency (MATLAB's default is -3). Inf when it never does.
static bool bandwidth(const ICoreTransferFunction& model, const double& dropDb,
double& frequency, std::string* whyNot = nullptr);
// The classical stability margins (toolbox row T1.20), MATLAB's
// `[Gm, Pm, Wcg, Wcp] = margin(sys)` -- and the two frequencies are NOT
// the two margins' own: `Wcg` is where the PHASE crosses -180 (which is
// what the gain margin is read at) and `Wcp` is where the GAIN crosses 1.
// MATLAB names them for the margin they belong to, not for the crossing
// they are, and a reading that swaps them is self-consistent and wrong.
//
// Gm absolute, NOT dB -- the plot shows dB and the number does not
// Pm degrees, in (-180, 180]
//
// Inf/NaN when a crossing does not exist, which is MATLAB's answer too: a
// model whose phase never reaches -180 has an infinite gain margin and no
// frequency to report it at.
//
// ⚠ THESE ARE EXACT AND MATLAB'S ARE NOT, the same split T1.21 measured
// for the peak gain. A crossing is a polynomial root here: for a
// continuous model `Im(N(jw) conj(D(jw))) = w (B_N A_D - A_N B_D)` and
// `|N|^2 - |D|^2` are both real polynomials in u = w^2, and for a discrete
// one the same two conditions are polynomials in z through the reversed
// (reciprocal) polynomials. MATLAB searches instead, and `allmargin.m`
// says how hard: `rtol = 1e-3; % relative accuracy on computed
// crossings/margins`. Measured on R2026a: `margin(tf(1, [1 2 3 1]))`
// answers Gm = 5.0000187267445027 at w = 1.7320535105357167 where the
// exact answers are 5 and sqrt(3) = 1.7320508075688772 -- an error of
// 1.9e-5, from its search and not from this one. The parity rows carry a
// band that covers it and the row note records the numbers.
static bool stabilityMargins(const ICoreTransferFunction& model, double& gainMargin,
double& phaseCrossFrequency, double& phaseMargin,
double& gainCrossFrequency, std::string* whyNot = nullptr);
};
};ICorePidController.h#
src/ICoreBlocks/ICoreMath/ControlSystems/StaticMethods/ICorePidController.h
ICorePidController#
ICorePidController.h:41 · class · nested Target, Tuned · 0 declaration(s)
MATLAB's pid and pidstd as TRANSFER FUNCTIONS -- toolbox row T1.30.
class ICorePidController {
public:
// Kp + Ki/s + Kd*s/(Tf*s + 1), or its ForwardEuler discretization when
// `Ts` is positive. `num` and `den` come back in MATLAB's descending
// order, normalised so the denominator is monic -- which is what MATLAB's
// own `tf(C)` answers.
//
// A negative `Tf`, or a `Ts` that is zero, is refused: MATLAB refuses
// both by name.
static bool parallelForm(const double& kp, const double& ki, const double& kd,
const double& tf, const double& sampleTime,
std::vector<double>& num, std::vector<double>& den,
std::string* whyNot = nullptr);
// `pidstd`'s gains as `pid`'s: Ki = Kp/Ti, Kd = Kp*Td, Tf = Td/N. An
// infinite Ti is "no integrator" (Ki = 0) and an infinite N is "no
// derivative filter" (Tf = 0), which is how MATLAB spells both.
static bool standardToParallel(const double& kp, const double& ti, const double& td,
const double& n, double& ki, double& kd, double& filter,
std::string* whyNot = nullptr);
// `pid(sys)` -- read the four gains back out of a transfer function that
// IS a parallel-form PID, and refuse one that is not.
//
// ⚠ THE REFUSAL IS THE POINT OF THIS FUNCTION. Every second-order model
// over `s` has SOME Kp/Ki/Kd that reproduces it, so a conversion that
// just solved the three equations would turn any plant into a controller.
// What makes a model a PID is its DENOMINATOR: `s` for the unfiltered
// form and `Tf*s^2 + s` for the filtered one, and nothing else.
static bool fromTransferFunction(const std::vector<double>& num,
const std::vector<double>& den,
const double& sampleTime, double& kp, double& ki,
double& kd, double& filter,
std::string* whyNot = nullptr);
// ---- pidtune (T1.31) --------------------------------------------------
//
// ⚠ THIS IS THE ONE ROW ON ITS BOARD THAT IS A DIFFERENT ALGORITHM BY
// DESIGN, and the row says so. MATLAB's `pidtune` is a proprietary
// loop-shaping search: `pidtune_` is not a `.m` in the installation, and
// what IS published is the CONTRACT -- "PIDTUNE tries to enforce a phase
// margin greater or equal to this value", default 60 degrees
// (`pidtuneOptions.m`). So the gains this answers are not MATLAB's and are
// not claimed to be; what both sides must agree on is the property, which
// is what T0.10 and X4 exist for.
//
// THE RULE, in full, because a documented rule is the whole difference
// between a design and a guess:
//
// 1. At the crossover `wc`, the loop L = C*G must have |L| = 1 and
// angle(L) = -180 + PM. Those two conditions fix the controller's
// response there completely: C(jwc) = (1/|G|) * exp(j*theta) with
// theta = -180 + PM - angle(G).
// 2. A structure can only supply theta from a RANGE, and the range is
// the structure: P supplies 0 exactly, I supplies angle(1/(jw))
// exactly, PI spans (-90, 0), PD spans (0, 90) and PID spans
// (-90, 90). A theta outside its structure's range is refused BY
// NAME, with the range and the margin that IS reachable -- MATLAB
// would silently do its best and answer a controller that misses the
// target it was given.
// 3. Inside the range the gains follow with no freedom left, except for
// PID, which has three gains against two conditions. The free
// parameter is spent on the CONTROLLER'S TWO ZEROS, which are placed
// COINCIDENT: C = Kd*(s + wc/b)^2/s, whose phase at jwc is
// 2*atan(b) - 90 and whose magnitude is Kd*wc*(1 + 1/b^2). So `b` is
// exactly the phase the design needs and nothing is left to choose.
// 4. When the caller names no `wc`, it is the frequency at which the
// required theta sits in the MIDDLE of the structure's range -- the
// most robust place to put it, and a rule rather than a search.
//
// ⚠ A DISCRETE PLANT IS NOT THE CONTINUOUS DESIGN WITH s REPLACED, for the
// same reason T1.30's discrete PID is not: MATLAB's default IFormula and
// DFormula are ForwardEuler, so the integrator is Ts/(z-1) and the
// derivative (z-1)/Ts, and neither has the phase 1/(jw) and jw have. Both
// domains run through ONE path here -- the two basis responses are
// evaluated numerically at wc and everything above is written over them --
// which is why step 3's closed form is solved by a bisection on `b` rather
// than by the formula the continuous case would allow.
enum class Form { P, I, PI, PD, PID };
// `crossover` <= 0 asks for step 4's automatic choice.
struct Target {
double phaseMarginDegrees = 60.0;
double crossover = 0.0;
};
// The gains, and the two numbers `[C, info] = pidtune(...)` reports beside
// them. ⚠ `phaseMargin` is the margin the design PUTS AT `crossover`: the
// target when the crossover was chosen here, and whatever the clamp left
// when the caller named one. It is not read back off the finished loop --
// `margin(C*G)` is what does that, it is exact on this console (T1.20),
// and the two agree whenever the loop crosses gain 1 only once.
struct Tuned {
double kp = 0.0;
double ki = 0.0;
double kd = 0.0;
double filter = 0.0; // always 0: 'pidf' is out of this row's scope
double crossover = 0.0;
double phaseMargin = 0.0;
};
static bool tune(const ICoreTransferFunction& plant, const Form& form, const Target& target,
Tuned& out, std::string* whyNot = nullptr);
};
};ICorePolePlacement.h#
src/ICoreBlocks/ICoreMath/ControlSystems/StaticMethods/ICorePolePlacement.h
ICorePolePlacement#
ICorePolePlacement.h:35 · class · 0 declaration(s)
State-feedback pole placement (the toolbox board's T1.25): the gain K that makes eig(A - B*K) the spectrum a caller asks for.
class ICorePolePlacement {
public:
// acker(A, B, p). B must have exactly one column. Answers a 1 x n row.
static bool ackermann(const ICoreMatrix& A, const ICoreMatrix& B,
const std::vector<std::complex<double>>& poles,
ICoreMatrix& gainOut, std::string* whyNot);
// place(A, B, p). Answers an (inputs) x n gain.
//
// `precisionOut`, when given, receives MATLAB's PREC -- the number of
// accurate decimal digits in the closed-loop poles, min(15, -log10(rel))
// over the non-zero requested poles, which is what MATLAB's second output
// is and what its 10% warning is measured from.
//
// ⚠ One step is not MATLAB's: where place.m maintains the QR of the
// eigenvector basis with rank-one qrinsert/qrdelete updates, this
// RECOMPUTES the factorization of the reduced basis each time. The
// quantity either route produces is the same one -- a unit vector
// orthogonal to every column but the one being improved -- and it enters
// the answer only through a subspace projection whose sign and scale are
// divided out, so the assigned spectrum is unaffected. It is recorded
// because it is the reason a MULTI-INPUT gain may differ from MATLAB's in
// the last digits where a single-input one (which never enters the
// refinement loop) does not.
static bool robust(const ICoreMatrix& A, const ICoreMatrix& B,
const std::vector<std::complex<double>>& poles,
ICoreMatrix& gainOut, double* precisionOut, std::string* whyNot);
};
};ICoreRiccatiDesign.h#
src/ICoreBlocks/ICoreMath/ControlSystems/StaticMethods/ICoreRiccatiDesign.h
ICoreRiccatiDesign#
ICoreRiccatiDesign.h:28 · class · 0 declaration(s)
The three answers MATLAB's Riccati solvers give (the toolbox board's T1.28): the stabilizing solution X, the state-feedback GAIN it implies, and the CLOSED-LOOP eigenvalues that gain produces.
class ICoreRiccatiDesign {
public:
// The continuous equation A'X + XA - XBR^-1B'X + Q = 0, and what MATLAB
// reads off it:
// K = R^-1 (B'X) the state-feedback gain
// L = eig(A - B*K) the closed-loop eigenvalues, a COLUMN
//
// R is the identity when a caller passes none -- MATLAB's own default for
// the three-argument care(A, B, Q), verified against R2026a.
//
// Answers false with `whyNot` set when the arguments do not conform or the
// solver could not find a stabilizing solution; the outputs are then
// untouched.
//
// `cross` is MATLAB's N, the term that couples the state and the input in
// the cost (2x'Nu). It is optional -- a null pointer is N = 0, which is
// what every three- and four-argument spelling means -- and it is folded
// in by the standard reduction rather than by a second solver:
//
// A - B R^-1 N' in place of A, Q - N R^-1 N' in place of Q
//
// so the solution X is the SAME solution and only the gain carries the
// extra N': K = R^-1 (B'X + N'). Doing it any other way would give the
// family two Riccati solvers to disagree with each other.
static bool continuousSolution(const ICoreMatrix& A, const ICoreMatrix& B,
const ICoreMatrix& Q, const ICoreMatrix& R,
ICoreMatrix& solutionOut, ICoreMatrix& gainOut,
ICoreComplexMatrix& closedLoopOut, std::string* whyNot,
const ICoreMatrix* cross = nullptr);
// The discrete equation A'XA - X - A'XB(R + B'XB)^-1B'XA + Q = 0, whose
// gain carries the extra B'XB term the continuous one does not:
// K = (B'XB + R)^-1 (B'XA)
// L = eig(A - B*K)
static bool discreteSolution(const ICoreMatrix& A, const ICoreMatrix& B,
const ICoreMatrix& Q, const ICoreMatrix& R,
ICoreMatrix& solutionOut, ICoreMatrix& gainOut,
ICoreComplexMatrix& closedLoopOut, std::string* whyNot,
const ICoreMatrix* cross = nullptr);
// ---- the estimator half (the toolbox board's T1.27) ------------------
//
// The Kalman gain is the LQR gain of the DUAL system, and that is not an
// analogy: `kalman` and `lqe` solve the very same Riccati equation with
// (A', C') in place of (A, B), and MATLAB's own kalman.m calls
// `icare(A', C', ...)` to get it. So this is the two solvers above read
// once through the transpose rather than a third solver -- which is what
// keeps a filter designed by `kalman` and a regulator designed by `lqr`
// from disagreeing about one equation.
//
// continuous A P + P A' - (P C' + N)(R)^-1(C P + N') + Q = 0
// L = (P C' + N) R^-1
// discrete P = A P A' - (A P C' + N)(C P C' + R)^-1(C P A' + N') + Q
// L = (A P C' + N)(C P C' + R)^-1
//
// Q, R and N here are the AGGREGATE covariances a plant's noise inputs
// imply (kalman's Qbar = G Qn G' and so on), not the user's Qn/Rn: forming
// them is the caller's job because only the caller knows which of the
// plant's inputs are noise.
//
// `covarianceOut` is P, `gainOut` is L (an Nx x Ny matrix, NOT its
// transpose -- MATLAB reports the gain as it multiplies the innovation),
// and `spectrumOut` is eig(A - L C), the estimator's own poles.
static bool estimatorSolution(const ICoreMatrix& A, const ICoreMatrix& C,
const ICoreMatrix& Q, const ICoreMatrix& R,
const bool& discrete, ICoreMatrix& covarianceOut,
ICoreMatrix& gainOut, ICoreComplexMatrix& spectrumOut,
std::string* whyNot, const ICoreMatrix* cross = nullptr);
};
};ICoreRootLocus.h#
src/ICoreBlocks/ICoreMath/ControlSystems/StaticMethods/ICoreRootLocus.h
ICoreRootLocus#
ICoreRootLocus.h:8 · class · 1 declaration(s)
Classical root locus under unity negative feedback.
class ICoreRootLocus {
public:
// For each gain K in `gains` (a vector), computes the closed-loop poles
// of the unity-negative-feedback system K*openLoopTf / (1 + K*openLoopTf)
// -- i.e. the roots of den(openLoopTf) + K*num(openLoopTf) -- and returns
// them all stacked as an (numGains * order) x 2 [real, imag] matrix.
// Poles are not ordered/matched into continuous branches across gains
// (a standard simplification for a first-pass locus plotter).
static ICoreMatrix computePoles(const ICoreTransferFunction& openLoopTf, const ICoreMatrix& gains);
};
};ICoreSS2TF.h#
src/ICoreBlocks/ICoreMath/ControlSystems/StaticMethods/ICoreSS2TF.h
ICoreSS2TF#
ICoreSS2TF.h:9 · class · 1 declaration(s)
class ICoreSS2TF {
public:
static ICoreArray<ICoreTransferFunction> ss2tf(const ICoreStateSpace& stateSpace);
};
};ICoreStateSpaceDiscretization.h#
src/ICoreBlocks/ICoreMath/ControlSystems/StaticMethods/ICoreStateSpaceDiscretization.h
The FOH discretization's STATE SHIFT, G2 = integral_0^Ts e^{A(Ts-s)} B s ds / Ts.
A first-order hold has no causal realization as it stands -- x[k+1] depends on u[k+1] -- so discretize_FOH() shifts the state by this matrix (xbar = x - G2 u) to get one, which is MATLAB's own move (
utDiscretizeFOH.m: "xd[k] = xc(k*Ts) + G * [u(...)]"). The shift is invisible in y, because D absorbs it, and it is exactly what a caller simulating withlsimneeds to hand back the CONTINUOUS state trajectory: xc[k] = xbar[k] + G2 u[k], and the initial condition goes the other way, xbar[0] = x0 - G2 u[0] (T1.15). Answered here rather than recomputed by the caller so that one augmented matrix exponential defines both halves.
ICoreStateSpaceDiscretization#
ICoreStateSpaceDiscretization.h:8 · class · 2 declaration(s)
class ICoreStateSpaceDiscretization {
public:
static ICoreStateSpace discretizeStateSpace(const ICoreStateSpace& originalContStateSpace, const double& Ts, const std::string& method = "ZOH");
static ICoreStateSpace deDiscretizeStateSpace(const ICoreStateSpace &originalDiscStateSpace, const std::string &method = "ZOH");
// The FOH discretization's STATE SHIFT, G2 = integral_0^Ts e^{A(Ts-s)} B s ds / Ts.
//
// A first-order hold has no causal realization as it stands -- x[k+1]
// depends on u[k+1] -- so discretize_FOH() shifts the state by this
// matrix (xbar = x - G2 u) to get one, which is MATLAB's own move
// (`utDiscretizeFOH.m`: "xd[k] = xc(k*Ts) + G * [u(...)]"). The shift is
// invisible in y, because D absorbs it, and it is exactly what a caller
// simulating with `lsim` needs to hand back the CONTINUOUS state
// trajectory: xc[k] = xbar[k] + G2 u[k], and the initial condition goes
// the other way, xbar[0] = x0 - G2 u[0] (T1.15). Answered here rather
// than recomputed by the caller so that one augmented matrix exponential
// defines both halves.
static ICoreMatrix firstOrderHoldStateShift(const ICoreStateSpace& continuousStateSpace,
const double& Ts);
};
};ICoreStructuralForms.h#
src/ICoreBlocks/ICoreMath/ControlSystems/StaticMethods/ICoreStructuralForms.h
ICoreStructuralForms#
ICoreStructuralForms.h:18 · class · 0 declaration(s)
The structural questions about a realization (the toolbox board's T1.24): which states an input can reach, which states an output can see, how far each one is driven or observed, and the change of ...
class ICoreStructuralForms {
public:
// ctrb(A, B) = [B, A*B, ..., A^(n-1)*B] and obsv(A, C) = [C; C*A; ...;
// C*A^(n-1)]. Built by repeated multiplication rather than by powers of A:
// the same matrix, one multiplication per block instead of one power.
static ICoreMatrix controllability(const ICoreMatrix& A, const ICoreMatrix& B,
std::string* whyNot);
static ICoreMatrix observability(const ICoreMatrix& A, const ICoreMatrix& C,
std::string* whyNot);
// The controllability staircase: an ORTHOGONAL T with
// Abar = T*A*T', Bbar = T*B, Cbar = C*T', in which the uncontrollable
// modes sit in the leading block and the controllable ones in the trailing
// one. `blocksOut` is MATLAB's `k`, the number of controllable states
// found at each iteration -- sum(k) is the controllable rank, which is the
// number this form exists to make readable.
//
// `tolerance` at or below zero asks for MATLAB's own default,
// n * norm(A, 1) * eps.
static bool staircase(const ICoreMatrix& A, const ICoreMatrix& B, const ICoreMatrix& C,
const double& tolerance, ICoreMatrix& aBarOut, ICoreMatrix& bBarOut,
ICoreMatrix& cBarOut, ICoreMatrix& transformOut,
ICoreMatrix& blocksOut, std::string* whyNot);
// The observability staircase, MATLAB's obsvf: the controllability
// staircase of the DUAL (A', C', B'), transposed back.
static bool dualStaircase(const ICoreMatrix& A, const ICoreMatrix& B, const ICoreMatrix& C,
const double& tolerance, ICoreMatrix& aBarOut,
ICoreMatrix& bBarOut, ICoreMatrix& cBarOut,
ICoreMatrix& transformOut, ICoreMatrix& blocksOut,
std::string* whyNot);
// gram(sys, 'c') and gram(sys, 'o') -- the energy each mode takes to reach
// and the energy each one produces. They are Lyapunov solutions and
// nothing more, which is why they are three lines here rather than a
// solver of their own:
//
// continuous A*Wc + Wc*A' + B*B' = 0 A'*Wo + Wo*A + C'*C = 0
// discrete A*Wc*A' - Wc + B*B' = 0 A'*Wo*A - Wo + C'*C = 0
//
// A gramian exists only for a STABLE A (MATLAB says so and refuses), and
// the Lyapunov solver's own singular-operator report is what says it here.
static ICoreMatrix gramian(const ICoreMatrix& A, const ICoreMatrix& BorC,
const bool& controllability, const bool& discrete,
std::string* whyNot);
// ---- the realization forms (the toolbox board's T1.29) ---------------
//
// Four ways to write the SAME system down. Every one of them is a change
// of basis on (A, B, C) -- so the transfer function is untouched by all
// four -- and what differs is which question the new basis makes readable:
// an arbitrary one (ss2ss), the characteristic polynomial (companion), the
// energy each mode carries (balanced), and then a smaller model built by
// dropping states the balanced form says carry almost none (modred).
//
// MATLAB's own convention for the transformation is xbar = T x, so
// Abar = T A T^-1, Bbar = T B, Cbar = C T^-1
// and `ss2ss(sys, T)`, `canon`'s second output and `balreal`'s third all
// mean that same T. It is the one thing to get right: the transpose
// convention answers a plausible model for a symmetric T and a wrong one
// for every other.
// ss2ss(sys, T): the basis the caller names. Exact -- three products.
static bool changeOfBasis(const ICoreMatrix& A, const ICoreMatrix& B, const ICoreMatrix& C,
const ICoreMatrix& T, ICoreMatrix& aOut, ICoreMatrix& bOut,
ICoreMatrix& cOut, std::string* whyNot);
// canon(sys, 'companion'): the basis in which the last column of A is the
// characteristic polynomial and B is e1. `canon.m`'s own construction --
// T = ctrb(A, B(:, 1)) and A_c = T \ A * T -- and the answer is
// INVARIANT under a prior diagonal scaling of the state, which is why the
// scaling MATLAB applies before it (ltipack.xscale) does not have to be
// reproduced here to get MATLAB's digits.
//
// `transformOut` is MATLAB's second output, the INVERSE of the internal T
// (canon.m returns `inv(T)`, "to be compatible with ss2ss").
static bool companionForm(const ICoreMatrix& A, const ICoreMatrix& B, const ICoreMatrix& C,
ICoreMatrix& aOut, ICoreMatrix& bOut, ICoreMatrix& cOut,
ICoreMatrix& transformOut, std::string* whyNot);
// balreal(sys): the basis in which the two gramians are EQUAL and diagonal,
// and that diagonal -- the Hankel singular values -- is `hankelOut`.
// Square-root balancing, which is balreal.m's own:
//
// Wc = Rr' Rr, Wo = Ro' Ro (Cholesky factors)
// Ro Rr' = U S V' (one SVD)
// T = S^-1/2 U' Ro, T^-1 = Rr' V S^-1/2
//
// The singular values are invariant (they are sqrt(eig(Wc Wo)) in any
// basis), so `hankelOut` is comparable digit for digit; the BASIS is unique
// only up to the sign of each state, because U and V are.
//
// A plant must be STABLE (the gramians are Lyapunov solutions) and MINIMAL
// (a zero Hankel value has no square root to divide by). MATLAB answers for
// an unstable plant by splitting it first (stabsep) and reports Inf for the
// unstable modes; that split is not here, and an unstable plant is refused
// rather than balanced in a basis that means nothing.
static bool balancedRealization(const ICoreMatrix& A, const ICoreMatrix& B,
const ICoreMatrix& C, const bool& discrete,
ICoreMatrix& aOut, ICoreMatrix& bOut, ICoreMatrix& cOut,
ICoreMatrix& hankelOut, ICoreMatrix& transformOut,
ICoreMatrix& inverseOut, std::string* whyNot);
// modred(sys, elim, 'MatchDC' | 'Truncate'): the model that is left when
// the states in `eliminated` (zero-based here, MATLAB's `elim` is
// one-based) are removed. Two ways to remove them, and they are not two
// spellings of one answer:
//
// MatchDC xdot2 = 0 is solved for x2 and substituted, so the reduced
// model has the FULL model's DC gain exactly
// Truncate the rows and columns are simply deleted, and the DC gain
// moves
//
// In discrete time "DC" is z = 1 rather than s = 0, so the elimination runs
// on A22 - I; using A22 there answers a model whose DC gain is wrong and
// whose poles are plausible.
static bool eliminateStates(const ICoreMatrix& A, const ICoreMatrix& B, const ICoreMatrix& C,
const ICoreMatrix& D, const bool& discrete,
const std::vector<size_t>& eliminated, const bool& matchDC,
ICoreMatrix& aOut, ICoreMatrix& bOut, ICoreMatrix& cOut,
ICoreMatrix& dOut, std::string* whyNot);
};
};ICoreTimeResponse.h#
src/ICoreBlocks/ICoreMath/ControlSystems/StaticMethods/ICoreTimeResponse.h
ICoreTimeResponse#
ICoreTimeResponse.h:15 · class · 1 declaration(s)
Time-domain response of a transfer function to the standard test inputs.
class ICoreTimeResponse {
public:
enum class Input { Step, Impulse, Ramp };
// Returns an N x 2 [time, output] matrix, or an empty matrix with `error`
// set. `duration` is in seconds and must be > 0; `points` is the requested
// sample count (>= 2) and is what sets the step for a continuous tf. For a
// discrete tf the step is fixed by its own Ts, so `points` is ignored and
// the count follows from duration / Ts.
static ICoreMatrix compute(const ICoreTransferFunction& tf, const Input& input,
const double& duration, const size_t& points,
std::string& error);
// A duration to simulate over when the caller named none -- what makes a
// bare `step(G)` answer instead of erroring. Read off the poles: long
// enough for the slowest mode to settle, or to ring a few times when
// nothing decays. Always finite and > 0, for any tf including a default
// one, so a caller may pass it straight to compute().
static double suggestedDuration(const ICoreTransferFunction& tf);
static std::string inputName(const Input& input); // "Step", "Impulse", "Ramp"
};
};