 |
LINE Solver (C++)
Templated C++ port of the LINE queueing solver
|
|
LINE Solver (C++)
Templated C++ port of the LINE queueing solver
|
The min-normal closure as a DIFFERENTIAL-ALGEBRAIC system: solver_fluid_dae.m. More...
#include <algorithm>#include <cmath>#include <cstddef>#include <limits>#include <string>#include <vector>#include "line/solvers/fluid/fluid_moments.h"#include "line/solvers/fluid/fluid_nonhyperbolic.h"#include "line/solvers/fluid/solver_fluid.h"#include "line/util/error.h"#include "line/util/lstsq.h"#include "line/util/lu.h"#include "line/util/matrix.h"#include "rodas.hpp"Go to the source code of this file.
Classes | |
| struct | line::fluid::FluidDaeConstraints |
| Finite capacity regions as linear admission constraints on the fluid state. More... | |
| struct | line::fluid::FluidDaeGates |
| Which events each cap throttles, and what happens to the mass it stops. More... | |
| struct | line::fluid::FluidDaeStaging |
| The waiting room outside a capped region, as fluid coordinates. More... | |
| struct | line::fluid::FluidDaeLegs |
| The two legs of every event under the active caps, and the waiting rooms. More... | |
| struct | line::fluid::FluidDaeConservation |
| Population conservation, one row per CLOSED chain, in state space. More... | |
| struct | line::fluid::FluidDaeOptions |
| Options that only the DAE route reads. More... | |
| struct | line::fluid::FluidDaeSystem |
| Everything the residual needs, gathered so the Newton can stay generic. More... | |
| struct | line::fluid::FluidDaeNewtonInfo |
| What the Newton reports about where it stopped. More... | |
| struct | line::fluid::FluidDaeSwitch |
| What a hybrid transient did, beside the trajectory. More... | |
Namespaces | |
| namespace | line |
| namespace | line::fluid |
Functions | |
| double | line::fluid::fluid_dae_reach (const std::vector< double > &row, const std::vector< std::size_t > &coord_class, const std::vector< double > &njobs, std::size_t K) |
| The largest row x the population can produce, ignoring the coupling. | |
| template<typename T> | |
| FluidDaeConstraints | line::fluid::fluid_dae_constraints (const qn::NetworkStruct< T > &sn, const FluidMomentTerms &terms) |
| template<typename T> | |
| FluidDaeGates | line::fluid::fluid_dae_gates (const qn::NetworkStruct< T > &sn, const FluidMomentTerms &t, const FluidDaeConstraints &con) |
| FluidDaeStaging | line::fluid::fluid_dae_staging (const FluidMomentTerms &t, const FluidDaeConstraints &con) |
| void | line::fluid::fluid_dae_extend (FluidDaeConstraints &con, const FluidDaeStaging &stg, const FluidMomentTerms &t) |
| Extend every cap to the staging coordinates that hold mass INSIDE it. | |
| FluidDaeLegs | line::fluid::fluid_dae_legs (const FluidMomentTerms &t, const FluidDaeGates &gates, const FluidDaeStaging &stg, const FluidDaeConstraints &con, const std::vector< std::size_t > &active, const std::vector< double > &sg, const std::vector< double > &r, const std::vector< double > &mult, bool staged_flow) |
| Shared by the steady-state residual and the transient right-hand side so that the two solve the SAME model and not two spellings of it. | |
| template<typename T> | |
| FluidDaeConservation | line::fluid::fluid_dae_conservation (const qn::NetworkStruct< T > &sn, const FluidMomentTerms &terms, const FluidDaeStaging &stg) |
| The conserved chains as equations. | |
| FluidDaeOptions | line::fluid::fluid_dae_options (const FluidOptions &opt, const FluidDaeOptions &dopt) |
| The controls the DAE route actually reads: the struct a caller pinned, with whatever options.config set on top of it. | |
| std::vector< std::size_t > | line::fluid::fluid_dae_closable (const FluidMomentTerms &t) |
| Stations whose variance enters the drift. | |
| std::vector< double > | line::fluid::fluid_dae_sigma_from (const Matrix< double > &Sigma, const FluidMomentTerms &t, const std::vector< std::size_t > &cidx) |
| Project a state-level covariance onto the per-station variances the drift reads. | |
| Matrix< double > | line::fluid::fluid_dae_clamp_tangent (const FluidDaeConstraints &con, const std::vector< std::size_t > &active, const FluidMomentTerms &t) |
| Orthogonal projector onto the subspace the CLAMPING caps leave free: I - R'(RR')^-1 R over the covariance coordinates. | |
| bool | line::fluid::fluid_dae_residual (const FluidDaeSystem &sysd, const std::vector< double > &u, std::vector< double > &G, std::vector< double > *rates_out, bool rethrow=false, std::vector< double > *fire_out=nullptr) |
| The coupled algebraic system, stacked: drift, conservation, closure consistency. | |
| void | line::fluid::fluid_dae_project (std::vector< double > &u, std::size_t nfree) |
| Onto the feasible box: the state is free, the variances are not. | |
| FluidDaeNewtonInfo | line::fluid::fluid_dae_newton (const FluidDaeSystem &sysd, std::vector< double > &u, double tol, std::size_t maxit) |
| Damped PROJECTED Newton with a finite-difference Jacobian and an Armijo backtrack on the residual norm. | |
| std::vector< std::vector< double > > | line::fluid::fluid_dae_integrate (const FluidMomentTerms &terms, const FluidDaeConservation &cons, const FluidClosure &closure, const std::vector< double > &x0, const std::vector< double > &grid, double tol, bool withcov=false, const std::vector< std::size_t > &closable=std::vector< std::size_t >()) |
| Integrate the closure as an index-1 DAE, with RODAS, and report the state at every point of grid. | |
| std::vector< double > | line::fluid::fluid_dae_hold_multipliers (const FluidMomentTerms &terms, const FluidDaeGates &gates, const FluidDaeStaging &stg, const FluidDaeConstraints &con, const std::vector< std::size_t > &active, const std::vector< double > &x, const std::vector< double > &sg, const FluidClosure &cl, const std::vector< double > &m0, bool &ok) |
| The multipliers that hold the active caps at this state, by small Newton. | |
| std::vector< std::vector< double > > | line::fluid::fluid_dae_integrate_hybrid (const FluidMomentTerms &terms, const FluidDaeConservation &cons, const FluidDaeConstraints &con, const FluidDaeGates &gates, const FluidDaeStaging &stg, const FluidClosure &closure, const std::vector< double > &x0In, const std::vector< double > &grid, double tol, std::vector< FluidDaeSwitch > *switches=nullptr) |
| The transient UNDER CAPS: one index-1 DAE per segment, restarted at every located crossing. | |
| void | line::fluid::fluid_dae_metrics (const FluidMomentTerms &terms, const std::vector< double > &x, const FluidClosure &cl, Matrix< double > &QN, Matrix< double > &UN, Matrix< double > &RN, Matrix< double > &TN, const std::vector< double > *rates=nullptr) |
| The metrics of one state, read exactly as the steady-state table reads them. | |
| template<class T> | |
| FluidSolution | line::fluid::solver_fluid_dae (const qn::NetworkStruct< T > &sn, const FluidOptions &opt, const FluidDaeOptions &dopt_in=FluidDaeOptions()) |
| solver_fluid_dae.m: the min-normal closure solved as one system. | |
| template<class T> | |
| std::vector< FluidTranPoint > | line::fluid::solver_fluid_dae_transient (const qn::NetworkStruct< T > &sn, const FluidOptions &opt, double t_end, std::size_t points=101, const std::vector< double > &out_grid=std::vector< double >(), const FluidDaeOptions &dopt_in=FluidDaeOptions()) |
| @@SolverFLD/getTranAvg for the DAE route: the metrics ALONG the trajectory. | |
The min-normal closure as a DIFFERENTIAL-ALGEBRAIC system: solver_fluid_dae.m.
SOLVER_FLUID_MOMENTS already solves a differential system (the mean) coupled to an algebraic one (the covariance). It solves them by SUCCESSIVE SUBSTITUTION: integrate the mean to its fixed point at a held variance, solve the Lyapunov equation there, extract sigma2, repeat. This states the same closure as one system and solves it as one system. The closure itself is unchanged – the drift, the rate factors and the Lyapunov equation are taken from FLUID_MOMENT_TERMS and FLUID_LYAPUNOV untouched – only the way the coupled equations are discharged.
STEADY STATE is an algebraic system, not an integration:
0 = D r(x, sigma2) drift residual, nstate rows
0 = C x - Nchain population conservation, one row per
closed chain
0 = sigma2 - sigmaOf(x, sigma2) closure consistency, one row per
closable station
solved simultaneously by a damped projected Newton. Integrating a stable ODE until it stops moving is a poor way to solve f(x)=0: the cost is set by the slowest mode of the model rather than by the accuracy wanted, which is why the stiff models are expensive on the substitution route.
TRANSIENT is an index-1 DAE with a SINGULAR MASS MATRIX, integrated by the vendored RODAS (third_party/rodas.hpp):
d/dt x = D r(x, sigma2) differential rows 0 = C x - Nchain algebraic rows
RODAS IS THE DEFAULT AND THE ONLY INTEGRATOR HERE. The rest of the fluid solver runs LSODA (line/util/lsoda.h) so that MATLAB, the JAR, Python and C++ agree in the last digits, but LSODA integrates y' = f and cannot carry a singular mass matrix at all, so it is not a candidate for this path. RODAS is a Rosenbrock method of order (3)4 for M y' = f with singular M, and being a fixed sequence of six linear solves rather than an iteration it has no convergence history for a future port to diverge on – which is what makes it the right choice for a route that the other three codebases do not have yet.
WHY THE CONSTRAINT IS WRITTEN DIFFERENTLY IN THE TWO MODES. At a fixed point every flow already balances, so the differentiated form d/dt(Cx) = 0 is satisfied by anything and pins nothing; the steady state therefore uses C x = N directly. The transient needs the opposite: the mass matrix zeroes one row per chain and the constraint residual is written there, which is the index-1 form RODAS integrates.
WHAT THE ALGEBRAIC CONSTRAINT BUYS. Population conservation otherwise holds only to integrator tolerance: it is a consequence of the drift (the rows of D sum to zero on a closed chain), never an equation. Writing it as a constraint also makes the Newton system solvable, because the drift Jacobian is singular along exactly the conserved directions – the same singularity FLUID_LYAPUNOV works around by projecting onto range(D) – so the constraint rows supply the missing rank instead of a pseudo-inverse hiding it.
WHY THE COVARIANCE IS NOT A NEWTON UNKNOWN. Sigma is nstate^2 entries, so a Jacobian over it is quartic work – strictly worse than the cubic Lyapunov solves it would replace. Sigma is LINEAR in itself for a held x, so it is eliminated by one Lyapunov solve per residual evaluation and only sigma2, M numbers, joins x in the unknown vector.
Definition in file fluid_dae.h.
1.18.0