diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml new file mode 100644 index 000000000..003f9654e --- /dev/null +++ b/.github/workflows/Documentation.yml @@ -0,0 +1,24 @@ +name: Documentation + +on: + push: + branches: + - main + tags: '*' + pull_request: + +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@latest + with: + version: '1' + - name: Install dependencies + run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' + - name: Build and deploy + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key + run: julia --project=docs/ docs/make.jl \ No newline at end of file diff --git a/Project.toml b/Project.toml index a63db53c2..00cc60ac6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelingToolkitStandardLibrary" uuid = "16a59e39-deab-5bd0-87e4-056b12336739" authors = ["Chris Rackauckas and Julia Computing"] -version = "0.1.1" +version = "1.2.0" [deps] IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" @@ -12,10 +12,10 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] IfElse = "0.1" -ModelingToolkit = "5.26, 6, 7" +ModelingToolkit = "8" OffsetArrays = "1" OrdinaryDiffEq = "5.56, 6" -Symbolics = "0.1, 1, 2, 3, 4" +Symbolics = "4.2.2" julia = "1.6" [extras] diff --git a/README.md b/README.md index 4c53c3413..964398136 100644 --- a/README.md +++ b/README.md @@ -1,389 +1,67 @@ # ModelingToolkitStandardLibrary.jl [![CI](https://github.com/SciML/ModelingToolkitStandardLibrary.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/SciML/ModelingToolkitStandardLibrary.jl/actions/workflows/CI.yml) +[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](http://mtkstdlib.sciml.ai/stable/) +[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](http://mtkstdlib.sciml.ai/dev/) -A standard library of components to model the world and beyond. +The ModelingToolkit Standard Library is standard library of components to model the world and beyond. -- [Electrical modeling](#electrical-modeling) - - [Basic electrical components](#basic-electrical-components) - - [Ground](#ground) - - [Resistor](#resistor) - - [Capacitor](#capacitor) - - [Ideal electrical elements](#ideal-electrical-elements) - - [Short](#short) - - [IdealOpAmp](#idealopamp) - - [Sensors](#sensors) - - [CurrentSensor](#currentsensor) - - [PotentialSensor](#potentialsensor) - - [VoltageSensor](#voltagesensor) - - [Voltage/current sources](#voltagecurrent-sources) - - [ConstantVoltage](#constantvoltage) - - [StepVoltage](#stepvoltage) - - [SineVoltage](#sinevoltage) -- [Thermal modeling](#thermal-modeling) - - [Basic thermal components](#basic-thermal-components) - - [Thermal ground](#thermal-ground) - - [Heat capacitor](#heat-capacitor) - - [Thermal conductor](#thermal-conductor) - - [Thermal resistor](#thermal-resistor) - - [Convective conductor](#convective-conductor) - - [Convective resistor](#convective-resistor) - - [Thermal radiation](#thermal-radiation) - - [Thermal collector](#thermal-collector) - - [Thermal sensors](#thermal-sensors) - - [Temperature sensor](#temperature-sensor) - - [Relative temperature sensor](#relative-temperature-sensor) - - [Heat flow sensor](#heat-flow-sensor) - - [Thermal sources](#thermal-sources) - - [Fixed heat flow](#fixed-heat-flow) - - [Fixed temperature](#fixed-temperature) +## Installation -## Electrical modeling +Assuming that you already have Julia correctly installed, it suffices to import +ModelingToolkitStandardLibrary.jl in the standard way: -Currently, ModelingToolkitStandardLibrary contains basic electrical components, ideal electrical elements, sensors, and voltage/current sources. +```julia +import Pkg; Pkg.add("ModelingToolkitStandardLibrary") +``` -### Basic electrical components +## Tutorials and Documentation -#### Ground +For information on using the package, +[see the stable documentation](https://mtkstdlib.sciml.ai/stable/). Use the +[in-development documentation](https://mtkstdlib.sciml.ai/dev/) for the version of +the documentation, which contains the unreleased features. -**Function**: `Ground(;name)` +## Libraries -**Description**: Ground node with the potential of zero and connector `g`. -Note that specifying the macro `@named sys = Ground()` is equivalent to setting -`sys = Ground(;name, sys)`. Either method will suffice, and there is no need to -type the name twice. The same principle applies to the other electrical components. +The following are the constituant libraries of the ModelingToolkit Standard Library. -**Connectors**: -- `g` +- [Basic Blocks](http://mtkstdlib.sciml.ai/dev/API/blocks/) +- [Electrical Components](http://mtkstdlib.sciml.ai/dev/API/electrical/) +- [Magnetic Components](http://mtkstdlib.sciml.ai/dev/API/magnetic/) +- [Thermal Components](http://mtkstdlib.sciml.ai/dev/API/thermal/) -### Resistor +## Example -**Function**: `Resistor(;name, R = 1.0)` +The following is the [RC Circuit Demonstration](http://mtkstdlib.sciml.ai/dev/tutorials/rc_circuit/): -**Observables**: -- `R`: resistance (negative, zero, positive) +```julia +using ModelingToolkit, OrdinaryDiffEq, Plots +using ModelingToolkitStandardLibrary.Electrical -**States** -- `v`: voltage +R = 1.0 +C = 1.0 +V = 1.0 +@variables t +@named resistor = Resistor(R=R) +@named capacitor = Capacitor(C=C) +@named source = ConstantVoltage(V=V) +@named ground = Ground() -**Connectors**: -- positive pin -- negative pin +rc_eqs = [ + connect(source.p, resistor.p) + connect(resistor.n, capacitor.p) + connect(capacitor.n, source.n, ground.g) + ] -### Capacitor +@named rc_model = ODESystem(rc_eqs, t, systems=[resistor, capacitor, source, ground]) +sys = structural_simplify(rc_model) +prob = ODAEProblem(sys, Pair[], (0, 10.0)) +sol = solve(prob, Tsit5()) +plot(sol, vars = [capacitor.v,resistor.i], + title = "RC Circuit Demonstration", + labels = ["Capacitor Voltage" "Resistor Current"]) +savefig("plot.png") +``` -**Function**: `Capacitor(;name, C = 1.0)` - -**Observables**: -- `C`: capacitance (zero or positive) - -**Connectors**: -- positive pin -- negative pin - ---- - -### Ideal electrical elements - -#### Short - -**Function**: `Short()` - -**Description**: Short cut branch. - -**Connectors**: -- positive pin -- negative pin - -#### IdealOpAmp - -**Function**: `IdealOpAmp(;name)` - -**Description**: The ideal operational amplifier. - -**States**: -- `v1`: voltage of the left port -- `v2`: voltage of the right port -- `i1`: current of the left port -- `i2`: current of the right port - -**Connectors**: -- positive pin (left port) -- negative pin (left port) -- positive pin (right port) -- negative pin (right port) - ---- - -### Sensors - -#### CurrentSensor - -**Function**: `CurrentSensor(;name)` - -**States** -- `i`: current value from the positive to the negative pin - -**Connectors**: -- positive pin -- negative pin - - -#### PotentialSensor - -**Function**: `PotentialSensor(;name)` - -**Connectors**: -- pin (which is to be measured) - -#### VoltageSensor - -**Function**: `VoltageSensor(;name)` - -**States**: -- `v`: value of voltage between the two pins - -**Connectors**: -- positive pin -- negative pin - ---- - -### Voltage/current sources - -#### ConstantVoltage - -**Function**: `ConstantVoltage(;name, V = 1.0)` - -**Description**: The source for an ideal constant voltage. - -**Observables**: -- `V`: value of constant voltage - -**Connectors**: -- positive pin -- negative pin - -#### StepVoltage - -**Function**: `StepVoltage(;name, offset = 0.0, starttime = 0.0, height = 0.0)` - -**Description**: Step voltage source. - -**Observables**: -- `offset`: voltage offset -- `starttime`: time offset -- `height`: height of the step - -**Connectors**: -- positive pin -- negative pin - -#### SineVoltage - -**Function**: `SineVoltage(;name, offset = 0.0, amplitude = 0.0, frequency = 0.0, starttime = 0.0, phase = 0.0)` - -**Description**: Sine voltage source. - -**Observables**: -- `offset`: voltage offset -- `amplitude`: amplitude of the sine wave -- `frequency`: frequency of the sine wave -- `starttime`: time offset -- `phase`: phase of the sine wave - -**Connectors**: -- positive pin -- negative pin - -## Thermal modeling - -ModelingToolkitStandardLibrary contains basic thermal components for modeling heat transfer and fluid heat flow. - -### Basic thermal components - -#### Thermal ground - -**Function**: `ThermalGround(;name)` - -**Description**: Thermal port for 1-dimensional heat transfer with temperature set to zero. -Note that specifying the macro `@named sys = ThermalGround()` is equivalent to setting -`sys = ThermalGround(;name, sys)`. Either method will suffice, and there is no need to -type the name twice. The same principle applies to the other thermal components. - -#### Heat capacitor - -**Function**: `HeatCapacitor(;name, C = 1.0)` - -**Observables**: -- `C`: heat capacity (zero or positive) - -**State**: -- `T`: temperature (in Kelvin) - -**Connectors**: -- heat port - -#### Thermal conductor - -**Function**: `ThermalConductor(;name, G = 1.0)` - -**Observables**: -- `G`: thermal conductance - -**States**: -- `Q_flow`: heat flow rate -- `T`: temperature (in Kelvin) - -**Connectors**: -- two heat ports - -#### Thermal resistor - -**Function**: `ThermalResistor(;name, R = 1.0)` - -**Description**: The model operates on the same principle as `ThermalConductor`, but relies on thermal resistance instead of thermal conductance. - -**Observables**: -- `R`: thermal resistance - -**States**: -- `Q_flow`: heat flow rate -- `T`: temperature (in Kelvin) - -**Connectors**: -- two heat ports - -#### Convective conductor - -**Function**: `ConvectiveConductor(;name; G = 1.0)` - -**Description**: Model of linear heat convection. - -**Observables**: -- `G`: convective thermal conductance - -**States**: -- `Q_flow`: heat flow rate -- `dT`: temperature difference (in Kelvin) - -**Connectors**: -- two heat ports (for modeling of the fluid flow over the solid) - -#### Convective resistor - -**Function**: `ConvectiveResistor(;name; R = 1.0)` - -**Description**: Model of linear heat convection. Works analogously to the above model, but relies on convective thermal resistance instead of convective thermal conductance. - -**Observables**: -- `R`: convective thermal resistance - -**States**: -- `Q_flow`: heat flow rate -- `dT`: temperature difference (in Kelvin) - -**Connectors**: -- two heat ports (for modeling of the fluid flow over the solid) - -#### Thermal radiation - -**Function**: `BodyRadiation(;name; G = 1.0)` - -**Description**: Thermal radiation model. - -**Observables**: -- `G`: net radiation conductance between two surfaces -- Stefan-Boltzmann constant - -**States**: -- `Q_flow`: heat flow rate - -**Connectors**: -- two heat ports - -#### Thermal collector - -**Function**: `ThermalCollector(;name, N = 1)` - -**Description**: A model for collecting the heat flows from multiple heatports -to a singular heatport. - -**Observables**: -- `N`: the number of input heatports - -**States**: -- `Q_flow`: heat flow rate -- `T`: temperature (in Kelvin) - -**Connectors**: -- `hp1...hpN`: the respective heatports -- `collector_port`: the target collector heatport - ---- - -### Thermal sensors - -#### Temperature sensor - -**Function**: `TemperatureSensor(;name)` - -**Description**: Ideal absolute temperature sensor which outputs the temperature (in Kelvin) of the connected port. - -**States**: -- `T`: temperature (in Kelvin) - -**Connectors**: -- heat port - -#### Relative temperature sensor - -**Function**: `RelativeTemperatureSensor(;name)` - -**Description**: The output of the sensor is the relative temperature, i.e., the difference of the two ports, given in Kelvin. - -**States**: -- `T`: temperature (in Kelvin) - -**Connectors**: -- two heat ports - -#### Heat flow sensor - -**Function**: `HeatFlowSensor(;name)` - -**Description**: The model monitors the heat flow rate of the component. Its output is positive when the direction the heat flow is from the first port to the second one. - -**States**: -- `Q_flow`: heat flow rate - -**Connectors**: -- two heat ports - ---- - -### Thermal sources - -#### Fixed heat flow - -**Function**: `FixedHeatFlow(;name, Q_flow=1.0, T₀=293.15, α=0.0)` - -**Observables**: -- `Q_flow`: the constant amount of heat flow rate -- `T₀`: the reference temperature -- `α`: this parameter simulates temperature-dependent loss (if the specified value is not 0) - -**Connectors**: -- heat port - -#### Fixed temperature - -**Function**: `FixedTemperature(;name, T = 0.0)` - -**Description**: The model defines a fixed temperature (in Kelvin) at a given port. - -**Observables**: -- `T`: temperature (in Kelvin) - -**Connectors**: -- heat port +![](https://user-images.githubusercontent.com/1814174/164912983-c3f73628-0e19-4e42-b085-4f62ba6f23d1.png) \ No newline at end of file diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 000000000..1bc937b99 --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,5 @@ +[deps] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" + +[compat] +Documenter = "0.26, 0.27" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 000000000..4bf65ef98 --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,47 @@ +using Documenter, ModelingToolkitStandardLibrary +using ModelingToolkitStandardLibrary.Blocks +using ModelingToolkitStandardLibrary.Mechanical +using ModelingToolkitStandardLibrary.Mechanical.Rotational +using ModelingToolkitStandardLibrary.Magnetic +using ModelingToolkitStandardLibrary.Magnetic.FluxTubes +using ModelingToolkitStandardLibrary.Electrical +using ModelingToolkitStandardLibrary.Thermal + +makedocs( + sitename="ModelingToolkitStandardLibrary.jl", + authors="Julia Computing", + clean=true, + doctest=false, + modules=[ModelingToolkitStandardLibrary, + ModelingToolkitStandardLibrary.Blocks, + ModelingToolkitStandardLibrary.Mechanical, + ModelingToolkitStandardLibrary.Mechanical.Rotational, + ModelingToolkitStandardLibrary.Magnetic, + ModelingToolkitStandardLibrary.Magnetic.FluxTubes, + ModelingToolkitStandardLibrary.Electrical, + ModelingToolkitStandardLibrary.Thermal], + + format=Documenter.HTML(assets=["assets/favicon.ico"], + canonical="https://mtkstdlib.sciml.ai/stable/"), + + pages=[ + "ModelingToolkitStandardLibrary.jl: A Standard Library for ModelingToolkit" => "index.md", + + "Tutorials" => [ + "RC Circuit" => "tutorials/rc_circuit.md" + ], + + "API" => [ + "Basic Blocks" => "API/blocks.md", + "Electrical Components" => "API/electrical.md", + "Magnetic Components" => "API/magnetic.md", + "Mechanical Components" => "API/mechanical.md", + "Thermal Components" => "API/thermal.md" + ], + ] +) + +deploydocs( + repo="github.com/SciML/ModelingToolkitStandardLibrary.jl"; + push_preview=true +) diff --git a/docs/src/API/blocks.md b/docs/src/API/blocks.md new file mode 100644 index 000000000..dd3c99e27 --- /dev/null +++ b/docs/src/API/blocks.md @@ -0,0 +1,76 @@ +# ModelingToolkitStandardLibrary: Blocks +```@meta +CurrentModule = ModelingToolkitStandardLibrary.Blocks +``` + +## Utility Blocks + +```@docs +RealInput +RealOutput +SISO +``` + +## Math Blocks + +```@docs +Gain +MatrixGain +Sum +Feedback +Add +Add3 +Product +Division +StaticNonLinearity +Abs +Sign +Sqrt +Sin +Cos +Tan +Asin +Acos +Atan +Atan2 +Sinh +Cosh +Tanh +Exp +Log +Log10 +``` + +## Source Blocks + +```@docs +Constant +Sine +Cosine +ContinuousClock +Ramp +Step +ExpSine +``` + +## Nonlinear Blocks + +```@docs +Limiter +DeadZone +SlewRateLimiter +``` + +## Continuous Blocks + +```@docs +Integrator +Derivative +FirstOrder +SecondOrder +StateSpace +PI +LimPI +PID +LimPID +``` \ No newline at end of file diff --git a/docs/src/API/electrical.md b/docs/src/API/electrical.md new file mode 100644 index 000000000..b57f30995 --- /dev/null +++ b/docs/src/API/electrical.md @@ -0,0 +1,52 @@ +# ModelingToolkitStandardLibrary: Electrical Components +```@meta +CurrentModule = ModelingToolkitStandardLibrary.Electrical +``` + +## Electrical Utilities + +```@docs +Pin +OnePort +``` + +## Analog Components + +```@docs +Capacitor +Ground +Inductor +Resistor +IdealOpAmp +``` + +## Analog Sensors + +```@docs +CurrentSensor +PotentialSensor +VoltageSensor +PowerSensor +MultiSensor +``` + +## Analogue Sources + +```@docs +ConstantVoltage +SineVoltage +StepVoltage +RampVoltage +SquareVoltage +TriangularVoltage +CosineVoltage +ExpSineVoltage +ConstantCurrent +SineCurrent +StepCurrent +RampCurrent +SquareCurrent +TriangularCurrent +CosineCurrent +ExpSineCurrent +``` \ No newline at end of file diff --git a/docs/src/API/magnetic.md b/docs/src/API/magnetic.md new file mode 100644 index 000000000..c3c7b0f6e --- /dev/null +++ b/docs/src/API/magnetic.md @@ -0,0 +1,34 @@ +# ModelingToolkitStandardLibrary: Magnetic Components + +## Flux Tubes +```@meta +CurrentModule = ModelingToolkitStandardLibrary.Magnetic.FluxTubes +``` + +### Flux Tube Utilities + +```@docs +PositiveMagneticPort +NegativeMagneticPort +TwoPort +``` + +### Basic Flux Tube Blocks + +```@docs +Ground +Idle +Short +Crossing +ConstantPermeance +ConstantReluctance +EddyCurrent +ElectroMagneticConverter +``` + +### Flux Tube Sources + +```@docs +ConstantMagneticPotentialDifference +ConstantMagneticFlux +``` \ No newline at end of file diff --git a/docs/src/API/mechanical.md b/docs/src/API/mechanical.md new file mode 100644 index 000000000..9ed292457 --- /dev/null +++ b/docs/src/API/mechanical.md @@ -0,0 +1,33 @@ +# ModelingToolkit Standard Library: Mechanical Components + +## Rotational Components +```@meta +CurrentModule = ModelingToolkitStandardLibrary.Mechanical.Rotational +``` + +### Rotational Utils + +```@docs +Flange +Support +PartialCompliantWithRelativeStates +PartialElementaryOneFlangeAndSupport2 +PartialElementaryTwoFlangesAndSupport2 +PartialCompliant +``` + +### Rotational Core Components + +```@docs +Fixed +Inertia +Spring +Damper +IdealGear +``` + +### Rotational Sources + +```@docs +Torque +``` \ No newline at end of file diff --git a/docs/src/API/thermal.md b/docs/src/API/thermal.md new file mode 100644 index 000000000..ebf5a467f --- /dev/null +++ b/docs/src/API/thermal.md @@ -0,0 +1,38 @@ +# ModelingToolkitStandardLibrary: Thermal Components +```@meta +CurrentModule = ModelingToolkitStandardLibrary.Thermal +``` + +## Thermal Utilities + +```@docs +HeatPort +Element1D +``` + +## Thermal Components + +```@docs +BodyRadiation +ConvectiveConductor +ConvectiveResistor +HeatCapacitor +ThermalConductor +ThermalResistor +ThermalCollector +``` + +## Thermal Sensors + +```@docs +RelativeTemperatureSensor +HeatFlowSensor +TemperatureSensor +``` + +## Thermal Sources + +```@docs +FixedHeatFlow +FixedTemperature +``` \ No newline at end of file diff --git a/docs/src/assets/favicon.ico b/docs/src/assets/favicon.ico new file mode 100644 index 000000000..3c6bd4703 Binary files /dev/null and b/docs/src/assets/favicon.ico differ diff --git a/docs/src/assets/logo.png b/docs/src/assets/logo.png new file mode 100644 index 000000000..6f4c3e261 Binary files /dev/null and b/docs/src/assets/logo.png differ diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 000000000..d98a23444 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,27 @@ +# ModelingToolkitStandardLibrary.jl + +ModelingToolkitStandardLibrary.jl is a standard library for the +[ModelingToolkit](https://mtk.sciml.ai/dev/) acasual modeling system. + +## Installation + +Assuming that you already have Julia correctly installed, it suffices to import +ModelingToolkitStandardLibrary.jl in the standard way: + +```julia +import Pkg; Pkg.add("ModelingToolkitStandardLibrary") +``` + +## Tutorials + +- [RC Circuit Tutorial](http://mtkstdlib.sciml.ai/dev/tutorials/rc_circuit/) + +## Libraries + +The following are the constituant libraries of the ModelingToolkit Standard Library. + +- [Basic Blocks](http://mtkstdlib.sciml.ai/dev/API/blocks/) +- [Electrical Components](http://mtkstdlib.sciml.ai/dev/API/electrical/) +- [Magnetic Components](http://mtkstdlib.sciml.ai/dev/API/magnetic/) +- [Mechanical Components](http://mtkstdlib.sciml.ai/dev/API/mechanical/) +- [Thermal Components](http://mtkstdlib.sciml.ai/dev/API/thermal/) diff --git a/docs/src/tutorials/rc_circuit.md b/docs/src/tutorials/rc_circuit.md new file mode 100644 index 000000000..0b7112181 --- /dev/null +++ b/docs/src/tutorials/rc_circuit.md @@ -0,0 +1,38 @@ +# RC Circuit Model + +This tutorial is a simplified version of the [RC circuit tutorial in the +ModelingToolkit.jl documentation](https://mtk.sciml.ai/dev/tutorials/acausal_components/). +In that tutorial, the full RC circuit is built from scratch. Here, we will use the +components of the `Electrical` model in the ModelingToolkit Standard Library to simply +connect pre-made components and simulate the model. + +```julia +using ModelingToolkit, OrdinaryDiffEq, Plots +using ModelingToolkitStandardLibrary.Electrical + +R = 1.0 +C = 1.0 +V = 1.0 +@variables t +@named resistor = Resistor(R=R) +@named capacitor = Capacitor(C=C) +@named source = ConstantVoltage(V=V) +@named ground = Ground() + +rc_eqs = [ + connect(source.p, resistor.p) + connect(resistor.n, capacitor.p) + connect(capacitor.n, source.n, ground.g) + ] + +@named rc_model = ODESystem(rc_eqs, t, systems=[resistor, capacitor, source, ground]) +sys = structural_simplify(rc_model) +prob = ODAEProblem(sys, Pair[], (0, 10.0)) +sol = solve(prob, Tsit5()) +plot(sol, vars = [capacitor.v,resistor.i], + title = "RC Circuit Demonstration", + labels = ["Capacitor Voltage" "Resistor Current"]) +savefig("plot.png") +``` + +![](https://user-images.githubusercontent.com/1814174/164912983-c3f73628-0e19-4e42-b085-4f62ba6f23d1.png) \ No newline at end of file diff --git a/plot.png b/plot.png new file mode 100644 index 000000000..0e387352d Binary files /dev/null and b/plot.png differ diff --git a/src/Blocks/Blocks.jl b/src/Blocks/Blocks.jl index 36bb597f0..4c5d75b56 100644 --- a/src/Blocks/Blocks.jl +++ b/src/Blocks/Blocks.jl @@ -1,29 +1,29 @@ """ The module `Blocks` contains common input-output components, referred to as blocks. - -In general, input-output blocks follow the convention -``` - ┌───────────┐ - u │ ẋ=f(x,u) │ y -────►│ y=g(x,u) ├────► - │ │ - └───────────┘ -``` -where `u` are inputs, `x` are state variables and `y` are outputs. `x,u,y` are all implemented as `@variables` internally, `u` are marked as `[input=true]` and `y` are marked `[output=true]`. """ module Blocks using ModelingToolkit, Symbolics, IfElse, OrdinaryDiffEq +using IfElse: ifelse @parameters t -Dₜ = Differential(t) +D = Differential(t) + +export RealInput, RealOutput, SISO +include("utils.jl") -export Gain, Sum +export Gain, Sum, MatrixGain, Feedback, Add, Product, Division +export Abs, Sign, Sqrt, Sin, Cos, Tan, Asin, Acos, Atan, Atan2, Sinh, Cosh, Tanh, Exp +export Log, Log10 include("math.jl") -export Saturation, DeadZone +export Constant, Sine, Cosine, ContinuousClock, Ramp, Step, ExpSine +include("sources.jl") + +export Limiter, DeadZone, SlewRateLimiter include("nonlinear.jl") -export Constant, Integrator, Derivative, FirstOrder, SecondOrder, PID, StateSpace +export Integrator, Derivative, FirstOrder, SecondOrder, StateSpace +export PI, LimPI, PID, LimPID include("continuous.jl") end \ No newline at end of file diff --git a/src/Blocks/continuous.jl b/src/Blocks/continuous.jl index eb62ba2e6..efda074a8 100644 --- a/src/Blocks/continuous.jl +++ b/src/Blocks/continuous.jl @@ -1,35 +1,22 @@ -# TODO: remove initial values for all inputs once IO handling in MTK is in place """ - Constant(val; name) - -Outputs a constant value `val`. -""" -function Constant(val; name) - @variables y(t)=val [output=true] - @parameters val=val - eqs = [ - y ~ val - ] - ODESystem(eqs, t, name=name) -end - -""" - Integrator(; k=1, name) + Integrator(;name, k=1, x_start=0.0) Outputs `y = ∫k*u dt`, corresponding to the transfer function `1/s`. """ -function Integrator(; k=1, name) - @variables x(t)=0 u(t)=0 [input=true] y(t)=0 [output=true] - @parameters k=k +function Integrator(;name, k=1, x_start=0.0) + @named siso = SISO() + @unpack u, y = siso + sts = @variables x(t)=x_start + pars = @parameters k=k eqs = [ - Dₜ(x) ~ k*u + D(x) ~ k * u y ~ x ] - ODESystem(eqs, t, name=name) + extend(ODESystem(eqs, t, sts, pars; name=name), siso) end """ - Derivative(; k=1, T, name) + Derivative(; name, k=1, T, x_start=0.0) Outputs an approximate derivative of the input. The transfer function of this block is ``` @@ -42,19 +29,27 @@ T 2 ⎛ 1⎞ and a state-space realization is given by `ss(-1/T, 1/T, -k/T, k/T)` where `T` is the time constant of the filter. A smaller `T` leads to a more ideal approximation of the derivative. + +# Parameters: +- `k`: Gain +- `T`: [s] Time constants (T>0 required; T=0 is ideal derivative block) +- `x_start`: Initial value of state """ -function Derivative(; k=1, T, name) - @variables x(t)=0 u(t)=0 [input=true] y(t)=0 [output=true] - @parameters T=T k=k +function Derivative(; name, k=1, T, x_start=0.0) + T > 0 || throw(ArgumentError("Time constant `T` has to be strictly positive")) + @named siso = SISO() + @unpack u, y = siso + sts = @variables x(t)=x_start + pars = @parameters T=T k=k eqs = [ - Dₜ(x) ~ (u - x) / T - y ~ (k/T)*(u - x) + D(x) ~ (u - x) / T + y ~ (k / T) * (u - x) ] - ODESystem(eqs, t, name=name) + extend(ODESystem(eqs, t, sts, pars; name=name), siso) end """ - FirstOrder(; k=1, T, name) + FirstOrder(; name, k=1, T, x_start=0.0) A first-order filter with a single real pole in `s = -T` and gain `k`. The transfer function is given by `Y(s)/U(s) = ` @@ -63,19 +58,27 @@ is given by `Y(s)/U(s) = ` ─────── sT + 1 ``` + +# Parameters: +- `k`: Gain +- `T`: [s] Time constants (T>0 required) +- `x_start`: Initial value of state """ -function FirstOrder(; k=1, T, name) - @variables x(t)=0 u(t)=0 [input=true] y(t) [output=true] - @parameters T=T k=k +function FirstOrder(; name, k=1, T, x_start=0.0) + T > 0 || throw(ArgumentError("Time constant `T` has to be strictly positive")) + @named siso = SISO() + @unpack u, y = siso + sts = @variables x(t)=x_start + pars = @parameters T=T k=k eqs = [ - Dₜ(x) ~ (-x + k*u) / T + D(x) ~ (k*u - x) / T y ~ x ] - ODESystem(eqs, t, name=name) + extend(ODESystem(eqs, t, sts, pars; name=name), siso) end """ - SecondOrder(; k=1, w, d, name) + SecondOrder(; name, k=1, w, d, x_start=0.0, xd_start=0.0) A second-order filter with gain `k`, a bandwidth of `w` rad/s and relative damping `d`. The transfer function is given by `Y(s)/U(s) = ` @@ -86,16 +89,161 @@ s² + 2d*w*s + w^2 ``` Critical damping corresponds to `d=1`, which yields the fastest step response without overshoot, d < 1` results in an under-damped filter while `d > 1` results in an over-damped filter. `d = 1/√2` corresponds to a Butterworth filter of order 2 (maximally flat frequency response). + +# Parameters: +- `k`: Gain +- `w`: Angular frequency +- `d`: Damping +- `x_start`: Initial value of state (output) +- `xd_start`: Initial value of derivative of state (output) """ -function SecondOrder(; k=1, w, d, name) - @variables x(t)=0 xd(t)=0 u(t)=0 [input=true] y(t) [output=true] - @parameters k=k w=w d=d +function SecondOrder(; name, k=1, w, d, x_start=0.0, xd_start=0.0) + @named siso = SISO() + @unpack u, y = siso + sts = @variables x(t)=x_start xd(t)=xd_start + pars = @parameters k=k w=w d=d eqs = [ - Dₜ(x) ~ xd - Dₜ(xd) ~ w*(w*(k*u - x) - 2*d*xd) + D(x) ~ xd + D(xd) ~ w*(w*(k*u - x) - 2*d*xd) y ~ x ] - ODESystem(eqs, t, name=name) + extend(ODESystem(eqs, t, sts, pars; name=name), siso) +end + +""" + PI(;name, k=1, T, x_start=0.0) + +Textbook version of a PI-controller without actuator saturation and anti-windup measure. + +# Parameters: +- `k`: Gain +- `T`: [s] Integrator time constant (T>0 required) +- `x_start`: Initial value for the integrator +""" +function PI(;name, k=1, T, x_start=0.0) + T > 0 || throw(ArgumentError("Time constant `T` has to be strictly positive")) + @named err_input = RealInput() # control error + @named ctr_output = RealOutput() # control signal + @named gainPI = Gain(k) + @named addPI = Add() + @named int = Integrator(k=1/T, x_start=x_start) + sys = [err_input, ctr_output, gainPI, addPI, int] + eqs = [ + connect(err_input, addPI.input1), + connect(addPI.output, gainPI.input), + connect(gainPI.output, ctr_output), + connect(err_input, int.input), + connect(int.output, addPI.input2), + ] + ODESystem(eqs, t, [], []; name=name, systems=sys) +end + +""" + PID(;name, k=1, Ti=false, Td=false, Nd=10, xi_start=0, xd_start=0) + +Text-book version of a PID-controller without actuator saturation and anti-windup measure. + +# Parameters: +- `k`: Gain +- `Ti`: [s] Integrator time constant (Ti>0 required). If set to false no integral action is used. +- `Td`: [s] Derivative time constant (Td>0 required). If set to false no derivative action is used. +- `Nd`: [s] Time constant for the derivative approximation (Nd>0 required; Nd=0 is ideal derivative). +- `x_start`: Initial value for the integrator. +- `xd_start`: Initial value for the derivative state. +""" +function PID(;name, k=1, Ti=false, Td=false, Nd=10, xi_start=0, xd_start=0) + with_I = !isequal(Ti, false) + with_D = !isequal(Td, false) + @named err_input = RealInput() # control error + @named ctr_output = RealOutput() # control signal + !isequal(Ti, false) && (Ti ≥ 0 || throw(ArgumentError("Ti out of bounds, got $(Ti) but expected Ti ≥ 0"))) + !isequal(Td, false) && (Td ≥ 0 || throw(ArgumentError("Td out of bounds, got $(Td) but expected Td ≥ 0"))) + Nd > 0 || throw(ArgumentError("Nd out of bounds, got $(Nd) but expected Nd > 0")) + + @named gainPID = Gain(k) + @named addPID = Add3() + if with_I + @named int = Integrator(k=1/Ti, x_start=xi_start) + else + @named Izero = Constant(k=0) + end + if with_D + @named der = Derivative(k=1/Td, T=1/Nd, x_start=xd_start) + else + @named Dzero = Constant(k=0) + end + sys = [err_input, ctr_output, gainPID, addPID] + if with_I + push!(sys, int) + else + push!(sys, Izero) + end + if with_D + push!(sys, der) + else + push!(sys, Dzero) + end + eqs = [ + connect(err_input, addPID.input1), + connect(addPID.output, gainPID.input), + connect(gainPID.output, ctr_output) + ] + if with_I + push!(eqs, connect(err_input, int.input)) + push!(eqs, connect(int.output, addPID.input2)) + else + push!(eqs, connect(err_input, Izero.input)) + push!(eqs, connect(Izero.output, addPID.input2)) + end + if with_D + push!(eqs, connect(err_input, der.input)) + push!(eqs, connect(der.output, addPID.input3)) + else + push!(eqs, connect(err_input, Dzero.input)) + push!(eqs, connect(Dzero.output, addPID.input3)) + end + ODESystem(eqs, t, [], []; name=name, systems=sys) +end + +""" + LimPI(;name, k=1, T, u_max=1, u_min=-u_max, Ta) + +Text-book version of a PI-controller with actuator saturation and anti-windup measure. + +# Parameters: +- `k`: Gain +- `T`: [s] Integrator time constant (T>0 required) +- `Ta`: [s] Tracking time constant (Ta>0 required) +- `x_start`: Initial value for the integrator +""" +function LimPI(;name, k=1, T, u_max=1, u_min=-u_max, Ta, x_start=0.0) + Ta > 0 || throw(ArgumentError("Time constant `Ta` has to be strictly positive")) + T > 0 || throw(ArgumentError("Time constant `T` has to be strictly positive")) + u_max ≥ u_min || throw(ArgumentError("u_min must be smaller than u_max")) + @named err_input = RealInput() # control error + @named ctr_output = RealOutput() # control signal + @named gainPI = Gain(k) + @named addPI = Add() + @named addTrack = Add() + @named int = Integrator(k=1/T, x_start=x_start) + @named limiter = Limiter(y_max=u_max, y_min=u_min) + @named addSat = Add(k1=1, k2=-1) + @named gainTrack = Gain(1/Ta) + sys = [err_input, ctr_output, gainPI, addPI, int, addTrack, limiter, addSat, gainTrack] + eqs = [ + connect(err_input, addPI.input1), + connect(addPI.output, gainPI.input), + connect(gainPI.output, limiter.input), + connect(limiter.output, ctr_output), + connect(limiter.input, addSat.input2), + connect(limiter.output, addSat.input1), + connect(addSat.output, gainTrack.input), + connect(err_input, addTrack.input1), + connect(gainTrack.output, addTrack.input2), + connect(addTrack.output, int.input), + connect(int.output, addPI.input2), + ] + ODESystem(eqs, t, [], []; name=name, systems=sys) end """ @@ -124,59 +272,98 @@ where the transfer function for the derivative includes additional filtering, se - `wd`: Set-point weighting in the derivative part. - `Nd`: Derivative limit, limits the derivative gain to Nd/Td. Reasonable values are ∈ [8, 20]. A higher value gives a better approximation of an ideal derivative at the expense of higher noise amplification. - `Ni`: `Ni*Ti` controls the time constant `Tₜ` of anti-windup tracking. A common (default) choice is `Tₜ = √(Ti*Td)` which is realized by `Ni = √(Td / Ti)`. Anti-windup can be effectively turned off by setting `Ni = Inf`. -`gains`: If `gains = true`, `Ti` and `Td` will be interpreted as gains with a fundamental PID transfer function on parallel form `ki=Ti, kd=Td, k + ki/s + kd*s` -""" -function PID(; k, Ti=false, Td=false, wp=1, wd=1, - Ni = Ti == 0 ? Inf : √(max(Td / Ti, 1e-6)), - Nd = 12, - y_max = Inf, - y_min = y_max > 0 ? -y_max : -Inf, - gains = false, - name -) +` `gains`: If `gains = true`, `Ti` and `Td` will be interpreted as gains with a fundamental PID transfer function on parallel form `ki=Ti, kd=Td, k + ki/s + kd*s` +""" +function LimPID(; name, k=1, Ti=false, Td=false, wp=1, wd=1, + Ni= Ti == 0 ? Inf : √(max(Td / Ti, 1e-6)), + Nd=10, + u_max=Inf, + u_min=u_max > 0 ? -u_max : -Inf, + gains=false, + xi_start=0.0, + xd_start=0.0, + ) + with_I = !isequal(Ti, false) + with_D = !isequal(Td, false) if gains Ti = k / Ti Td = Td / k end 0 ≤ wp ≤ 1 || throw(ArgumentError("wp out of bounds, got $(wp) but expected wp ∈ [0, 1]")) 0 ≤ wd ≤ 1 || throw(ArgumentError("wd out of bounds, got $(wd) but expected wd ∈ [0, 1]")) - Ti ≥ 0 || throw(ArgumentError("Ti out of bounds, got $(Ti) but expected Ti ≥ 0")) - Td ≥ 0 || throw(ArgumentError("Td out of bounds, got $(Td) but expected Td ≥ 0")) - y_max ≥ y_min || throw(ArgumentError("y_min must be smaller than y_max")) - - @variables x(t)=0 u_r(t)=0 [input=true] u_y(t)=0 [input=true] y(t) [output=true] e(t)=0 ep(t)=0 ed(t)=0 ea(t)=0 - + !isequal(Ti, false) && (Ti ≥ 0 || throw(ArgumentError("Ti out of bounds, got $(Ti) but expected Ti ≥ 0"))) + !isequal(Td, false) && (Td ≥ 0 || throw(ArgumentError("Td out of bounds, got $(Td) but expected Td ≥ 0"))) + u_max ≥ u_min || throw(ArgumentError("u_min must be smaller than u_max")) + Nd > 0 || throw(ArgumentError("Nd out of bounds, got $(Nd) but expected Nd > 0")) - @named D = Derivative(k = Td, T = Td/Nd) # NOTE: consider T = max(Td/Nd, 100eps()), but currently errors since a symbolic variable appears in a boolean expression in `max`. - if isequal(Ti, false) - @named I = Gain(false) + @named reference = RealInput() + @named measurement = RealInput() + @named ctr_output = RealOutput() # control signal + @named addP = Add(k1=wp, k2=-1) + @named gainPID = Gain(k) + @named addPID = Add3() + if with_I + @named addI = Add3(k1=1, k2=-1, k3=1) + @named int = Integrator(k=1/Ti, x_start=xi_start) + @named limiter = Limiter(y_max=u_max, y_min=u_min) + @named addSat = Add(k1=1, k2=-1) + @named gainTrack = Gain(1/(k * Ni)) else - @named I = Integrator(k = 1/Ti) + @named Izero = Constant(k=0) end - @named sat = Saturation(; y_min, y_max) - derivative_action = Td > 0 - @parameters k=k Td=Td wp=wp wd=wd Ni=Ni Nd=Nd # TODO: move this line above the subsystem definitions when symbolic default values for parameters works. https://github.com/SciML/ModelingToolkit.jl/issues/1013 - # NOTE: Ti is not included as a parameter since we cannot support setting it to false after this constructor is called. Maybe Integrator can be tested with Ti = false setting k to 0 with IfElse? - + if with_D + @named der = Derivative(k=1/Td, T=1/Nd, x_start=xd_start) + @named addD = Add(k1=wd, k2=-1) + else + @named Dzero = Constant(k=0) + end + + sys = [reference, measurement, ctr_output, addP, gainPID, addPID] + if with_I + push!(sys, [addI, int, limiter, addSat, gainTrack]...) + else + push!(sys, Izero) + end + if with_D + push!(sys, [addD, der]...) + else + push!(sys, Dzero) + end + eqs = [ - e ~ u_r - u_y # Control error - ep ~ wp*u_r - u_y # Control error for proportional part with setpoint weight - ea ~ sat.y - sat.u # Actuator error due to saturation - I.u ~ e + 1/(k*Ni)*ea # Connect integrator block. The integrator integrates the control error and the anti-wind up tracking. Note the apparent tracking time constant 1/(k*Ni), since k appears after the integration and 1/Ti appears in the integrator block, the final tracking gain will be 1/(Ti*Ni) - sat.u ~ derivative_action ? k*(ep + I.y + D.y) : k*(ep + I.y) # unsaturated output = P + I + D - y ~ sat.y + connect(reference, addP.input1), + connect(measurement, addP.input2), + connect(addP.output, addPID.input1), + connect(addPID.output, gainPID.input), ] - systems = [I, sat] - if derivative_action - push!(eqs, ed ~ wd*u_r - u_y) - push!(eqs, D.u ~ ed) # Connect derivative block - push!(systems, D) + if with_I + push!(eqs, connect(reference, addI.input1)) + push!(eqs, connect(measurement, addI.input2)) + push!(eqs, connect(gainPID.output, limiter.input)) + push!(eqs, connect(limiter.output, ctr_output)) + push!(eqs, connect(limiter.input, addSat.input2)) + push!(eqs, connect(limiter.output, addSat.input1)) + push!(eqs, connect(addSat.output, gainTrack.input)) + push!(eqs, connect(gainTrack.output, addI.input3)) + push!(eqs, connect(addI.output, int.input)) + push!(eqs, connect(int.output, addPID.input3)) + else + push!(eqs, connect(Izero.output, addPID.input3)) + end + if with_D + push!(eqs, connect(reference, addD.input1)) + push!(eqs, connect(measurement, addD.input2)) + push!(eqs, connect(addD.output, der.input)) + push!(eqs, connect(der.output, addPID.input2)) + else + push!(eqs, connect(Dzero.output, addPID.input2)) end - ODESystem(eqs, t, name=name, systems=systems) + + ODESystem(eqs, t, [], []; name=name, systems=sys) end """ - StateSpace(A, B, C, D=0; x0=zeros(size(A,1)), name) + StateSpace(A, B, C, D=0; x_start=zeros(size(A,1)), name) A linear, time-invariant state-space system on the form. ``` @@ -185,28 +372,26 @@ y = Cx + Du ``` Transfer functions can also be simulated by converting them to a StateSpace form. """ -function StateSpace(A, B, C, D=0; x0=zeros(size(A,1)), name) - nx = size(A,1) - nu = size(B,2) - ny = size(C,1) - if nx == 0 - length(C) == length(B) == 0 || throw(ArgumentError("Dimension mismatch between A,B,C matrices")) - return Gain(D; name=name) - end +function StateSpace(;A, B, C, D=nothing, x_start=zeros(size(A,1)), name) + nx, nu, ny = size(A,1), size(B,2), size(C,1) + size(A,2) == nx || error("`A` has to be a square matrix.") + size(B,1) == nx || error("`B` has to be of dimension ($nx x $nu).") + size(C,2) == nx || error("`C` has to be of dimension ($ny x $nx).") if B isa AbstractVector B = reshape(B, length(B), 1) end - if D == 0 + if isnothing(D) || iszero(D) D = zeros(ny, nu) + else + size(D) == (ny,nu) || error("`D` has to be of dimension ($ny x $nu).") end - @variables x[1:nx](t)=x0 u[1:nu](t)=0 [input=true] y[1:ny](t)=C*x0 [output=true] - x = collect(x) # https://github.com/JuliaSymbolics/Symbolics.jl/issues/379 - u = collect(u) - y = collect(y) - # @parameters A=A B=B C=C D=D # This is buggy - eqs = [ - Dₜ.(x) .~ A*x .+ B*u - y .~ C*x .+ D*u + @named input = RealInput(nin=nu) + @named output = RealOutput(nout=ny) + @variables x[1:nx](t)=x_start + # pars = @parameters A=A B=B C=C D=D # This is buggy + eqs = [ # FIXME: if array equations work + [Differential(t)(x[i]) ~ sum(A[i,k] * x[k] for k in 1:nx) + sum(B[i,j] * input.u[j] for j in 1:nu) for i in 1:nx]..., # cannot use D here + [output.u[j] ~ sum(C[j,i] * x[i] for i in 1:nx) + sum(D[j,k] * input.u[k] for k in 1:nu) for j in 1:ny]..., ] - ODESystem(eqs, t, name=name) + compose(ODESystem(eqs, t, vcat(x...), [], name=name), [input, output]) end diff --git a/src/Blocks/math.jl b/src/Blocks/math.jl index d0078b49c..86c7d382c 100644 --- a/src/Blocks/math.jl +++ b/src/Blocks/math.jl @@ -1,49 +1,281 @@ -import Symbolics.scalarize - """ - Gain(k; name) + Gain(k=1; name) + +Output the product of a gain value with the input signal. -Outputs `y = k*u`. `k` can be a scalar or an array. +# Parameters: +- `k`: Scalar gain """ function Gain(k=1; name) - @variables u(t)=0 [input=true] y(t) [output=true] - @parameters k=k + @named siso = SISO() + @unpack u, y = siso + pars = @parameters k=k eqs = [ - y ~ k*u + y ~ k * u ] - ODESystem(eqs, t, name=name) + extend(ODESystem(eqs, t, [], pars; name=name), siso) end -function Gain(K::AbstractArray; name) - ny,nu = size(K, 1), size(K, 2) - @variables u[1:nu](t)=0 [input=true] y[1:ny](t)=0 [output=true] - eqs = y .~ K*u - ODESystem(eqs, t, name=name) +""" + MatrixGain(K::AbstractArray; name) + +Output the product of a gain matrix with the input signal vector. + +# Parameters: +- `K`: Matrix gain +""" +function MatrixGain(K::AbstractArray; name) + nout, nin = size(K) + @named input = RealInput(;nin=nin) + @named output = RealOutput(;nout=nout) + eqs = [ + output.u[i] ~ sum(K[i,j] * input.u[j] for j in 1:nin) for i in 1:nout # FIXME: if array equations work + ] + compose(ODESystem(eqs, t, [], []; name=name), [input, output]) end """ Sum(n::Int; name) - Sum(k::AbstractVector; name) -Creates a summing block that sums `n` inputs, `y = sum(u[i] for i ∈ 1:n)`. -A vector of summing coefficients `k` can also be provided, i.e., `y = sum(k[i]u[i] for i ∈ 1:n)`. -A block that subtracts one signal from another can thus be created by `@named sub = Sum([1, -1])`. +Output the sum of the elements of the input port vector. + +# Parameters: +- `n`: Input port dimension """ function Sum(n::Int; name) - @variables u[1:n](t)=0 [input=true] y(t)=0 [output=true] - eqs = [y ~ scalarize(sum(u))] - ODESystem(eqs, t, name=name) + @named input = RealInput(;nin=n) + @named output = RealOutput() + eqs = [ + output.u ~ sum(input.u) + ] + compose(ODESystem(eqs, t, [], []; name=name), [input, output]) end + +""" + Feedback(;name) -function Sum(k::AbstractVector; name) - n = length(k) - @variables u[1:n](t)=0 [input=true] y(t)=0 [output=true] - eqs = [y ~ scalarize(sum(k[i]*u[i] for i ∈ 1:n))] - ODESystem(eqs, t, name=name) +Output difference between reference input (input1) and feedback input (input2). +""" +function Feedback(;name) + @named input1 = RealInput() + @named input2 = RealInput() + @named output = RealOutput() + eqs= [ + output.u ~ input1.u - input2.u + ] + return compose(ODESystem(eqs, t, [], []; name=name), input1, input2, output) end -function Product(n::Int=2; name) - @variables u[1:n](t)=0 [input=true] y(t)=0 [output=true] - eqs = [y ~ scalarize(prod(u))] - ODESystem(eqs, t, name=name) -end \ No newline at end of file +""" + Add(;name, k1=1, k2=1) + +Output the sum of the two scalar inputs. + +# Parameters: +- `k1`: Gain for first input +- `k2`: Gain for second input +""" +function Add(;name, k1=1, k2=1) + @named input1 = RealInput() + @named input2 = RealInput() + @named output = RealOutput() + pars = @parameters begin + k1=k1 + k2=k2 + end + eqs= [ + output.u ~ k1 * input1.u + k2 * input2.u + ] + return compose(ODESystem(eqs, t, [], pars; name=name), input1, input2, output) +end + +""" + Add(;name, k1=1, k2=1,k3=1) + +Output the sum of the three scalar inputs. + +# Parameters: +- `k1`: Gain for first input +- `k2`: Gain for second input +- `k3`: Gain for third input +""" +function Add3(;name, k1=1, k2=1, k3=1) + @named input1 = RealInput() + @named input2 = RealInput() + @named input3 = RealInput() + @named output = RealOutput() + pars = @parameters begin + k1=k1 + k2=k2 + k3=k3 + end + eqs= [ + output.u ~ k1 * input1.u + k2 * input2.u + k3 * input3.u + ] + return compose(ODESystem(eqs, t, [], pars; name=name), input1, input2, input3, output) +end + +""" + Product(;name) + +Output product of the two inputs. +""" +function Product(;name) + @named input1 = RealInput() + @named input2 = RealInput() + @named output = RealOutput() + eqs= [ + output.u ~ input1.u * input2.u + ] + return compose(ODESystem(eqs, t, [], []; name=name), input1, input2, output) +end + +""" + Division(;name) + +Output first input divided by second input. +""" +function Division(;name) + @named input1 = RealInput() + @named input2 = RealInput() + @named output = RealOutput() + eqs= [ + output.u ~ input1.u / input2.u + ] + return compose(ODESystem(eqs, t, [], []; name=name), input1, input2, output) +end + + +""" + StaticNonLinearity(func ;name) + +Applies the given function to the input. + +If the given function is not composed of simple core methods (e.g. sin, abs, ...), it has to be registered via `@register_symbolic func(u)` +""" +function StaticNonLinearity(func; name) + @named siso = SISO() + @unpack u, y = siso + eqs = [y ~ func(u)] + extend(ODESystem(eqs, t, [], []; name=name), siso) +end + +""" + Abs(;name) + +Output the absolute value of the input. +""" +Abs(;name) = StaticNonLinearity(abs; name) + +""" + Sign(;name) + +Output the sign of the input +""" +Sign(;name) = StaticNonLinearity(sign; name) + +""" + Sqrt(;name) + +Output the square root of the input (input >= 0 required). +""" +Sqrt(;name) = StaticNonLinearity(sqrt; name) + +""" + Sin(;name) + +Output the sine of the input. +""" +Sin(;name) = StaticNonLinearity(sin; name) + +""" + Cos(;name) + +Output the cosine of the input. +""" +Cos(;name) = StaticNonLinearity(cos; name) + +""" + Tan(;name) + +Output the tangent of the input. +""" +Tan(;name) = StaticNonLinearity(tan; name) + +""" + Asin(;name) + +Output the arc sine of the input. +""" +Asin(;name) = StaticNonLinearity(asin; name) + +""" + Acos(;name) + +Output the arc cosine of the input. +""" +Acos(;name) = StaticNonLinearity(acos; name) + +""" + Atan(;name) + +Output the arc tangent of the input. +""" +Atan(;name) = StaticNonLinearity(atan; name) + +""" + Atan2(;name) + +Output the arc tangent of the input. +""" +function Atan2(;name) + @named input1 = RealInput() + @named input2 = RealInput() + @named output = RealOutput() + eqs = [ + output.u ~ atan(input1.u, input2.u) + ] + compose(ODESystem(eqs, t, [], []; name=name), [input1, input2, output]) +end + +""" + Sinh(;name) + +Output the hyperbolic sine of the input. +""" +Sinh(;name) = StaticNonLinearity(sinh; name) + +""" + Cosh(;name) + +Output the hyperbolic cosine of the input. +""" +Cosh(;name) = StaticNonLinearity(cosh; name) + +""" + Tanh(;name) + +Output the hyperbolic tangent of the input. +""" +Tanh(;name) = StaticNonLinearity(tanh; name) + +""" + Exp(;name) + +Output the exponential (base e) of the input. +""" +Exp(;name) = StaticNonLinearity(exp; name) + +""" + Log(;name) + +Output the natural (base e) logarithm of the input. +""" +Log(;name) = StaticNonLinearity(log; name) + +""" + Log10(;name) + +Output the base 10 logarithm of the input. +""" +Log10(;name) = StaticNonLinearity(log10; name) \ No newline at end of file diff --git a/src/Blocks/nonlinear.jl b/src/Blocks/nonlinear.jl index 333ee7172..2a9ac5c0e 100644 --- a/src/Blocks/nonlinear.jl +++ b/src/Blocks/nonlinear.jl @@ -1,31 +1,29 @@ -const ie = IfElse.ifelse +_clamp(u, u_min, u_max) = max(min(u, u_max), u_min) +_dead_zone(u, u_min, u_max) = ifelse(u > u_max, u - u_max, ifelse(u < u_min, u - u_min, 0)) """ - Saturation(; y_max, y_min=-y_max, name) +Limit the range of a signal. -The `Saturation` IO block limits the output between `y_min` and `y_max`, equivalent to -`y ~ clamp(u, y_min, y_max)`. - -Keywords: limiter, sat, actuator model +# Parameters: +- `y_max`: Maximum of output signal +- `y_min`: Minimum of output signal """ -function Saturation(; y_max, y_min=y_max > 0 ? -y_max : -Inf, name) - if !ModelingToolkit.isvariable(y_max) - y_max ≥ y_min || throw(ArgumentError("y_min must be smaller than y_max")) - end - @variables u(t)=0 [input=true] y(t)=0 [output=true] - @parameters y_max=y_max y_min=y_min +function Limiter(;name, y_max, y_min=y_max > 0 ? -y_max : -Inf) + y_max ≥ y_min || throw(ArgumentError("`y_min` must be smaller than `y_max`")) + @named siso = SISO() + @unpack u, y = siso + pars = @parameters y_max=y_max y_min=y_min eqs = [ - # The equation below is equivalent to y ~ clamp(u, y_min, y_max) - y ~ ie(u > y_max, y_max, ie( (y_min < u) & (u < y_max), u, y_min)) + y ~ _clamp(u, y_min, y_max) ] - ODESystem(eqs, t, name=name) + extend(ODESystem(eqs, t, [], pars; name=name), siso) end """ DeadZone(; u_max, u_min=-u_max, name) -A dead zone is a band within which the output is zero. -Outside of the dead zone, the output changes linearly starting from zero at the band edge. +The DeadZone block defines a region of zero output. +If the input is within uMin ... uMax, the output is zero. Outside of this zone, the output is a linear function of the input with a slope of 1. ``` y▲ │ / @@ -37,14 +35,37 @@ Outside of the dead zone, the output changes linearly starting from zero at the / │ ``` """ -function DeadZone(; u_max, u_min=-u_max, name) +function DeadZone(; name, u_max, u_min=-u_max) if !ModelingToolkit.isvariable(u_max) - u_max ≥ u_min || throw(ArgumentError("u_min must be smaller than u_max")) + u_max ≥ u_min || throw(ArgumentError("`u_min` must be smaller than `u_max`")) end - @variables u(t)=0 [input=true] y(t)=0 [output=true] - @parameters u_max=u_max u_min=u_min + @named siso = SISO() + @unpack u, y = siso + pars = @parameters u_max=u_max u_min=u_min + eqs = [ + y ~ _dead_zone(u, u_min, u_max) + ] + extend(ODESystem(eqs, t, [], pars; name=name), siso) +end + +""" + SlewRateLimiter(;name, rising=1, falling=-rising, Td=0.001, y_start=0.0) + +Limits the slew rate of a signal. + +# Parameters: +- `Rising`: Maximum rising slew rate +- `falling`: Maximum falling slew rate +- `Td`: Derivative time constant +""" +function SlewRateLimiter(;name, rising=1, falling=-rising, Td=0.001, y_start=0.0) + rising ≥ falling || throw(ArgumentError("`rising` must be smaller than `falling`")) + Td > 0 || throw(ArgumentError("Time constant `Td` must be strictly positive")) + @named siso = SISO(y_start=y_start) + @unpack u, y = siso + pars = @parameters rising=rising falling=falling eqs = [ - y ~ ie(u > u_max, u-u_max, ie( u < u_min, u-u_min, 0)) + D(y) ~ max(min((u-y) / Td, rising), falling) ] - ODESystem(eqs, t, name=name) + extend(ODESystem(eqs, t, [], pars; name=name), siso) end diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl new file mode 100644 index 000000000..796227569 --- /dev/null +++ b/src/Blocks/sources.jl @@ -0,0 +1,154 @@ +""" +Generate constant signal. + +# Parameters: +- `k`: Constant output value +""" +function Constant(;name, k=1) + @named output = RealOutput() + pars = @parameters k=k + eqs = [ + output.u ~ k + ] + compose(ODESystem(eqs, t, [], pars; name=name), [output]) +end + +""" +Generate sine signal. + +# Parameters: +- `frequency`: [Hz] Frequency of sine wave +- `amplitude`: Amplitude of sine wave +- `phase`: [rad] Phase of sine wave +- `offset`: Offset of output signal +- `start_time`: [s] Output `y = offset` for `t < start_time` +""" +function Sine(;name, + frequency, + amplitude=1, + phase=0, + offset=0, + start_time=0) + + @named output = RealOutput() + pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase + eqs = [ + output.u ~ offset + ifelse(t < start_time, 0, amplitude* sin(2*pi*frequency*(t - start_time) + phase)) + ] + compose(ODESystem(eqs, t, [], pars; name=name), [output]) +end + +""" +Generate cosine signal. + +# Parameters: +- `frequency`: [Hz] Frequency of sine wave +- `amplitude`: Amplitude of sine wave +- `phase`: [rad] Phase of sine wave +- `offset`: Offset of output signal +- `start_time`: [s] Output `y = offset` for `t < start_time` +""" +function Cosine(;name, + frequency, + amplitude=1, + phase=0, + offset=0, + start_time=0) + + @named output = RealOutput() + pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase + eqs = [ + output.u ~ offset + ifelse(t < start_time, 0, amplitude* cos(2*pi*frequency*(t - start_time) + phase)) + ] + compose(ODESystem(eqs, t, [], pars; name=name), [output]) +end + +""" +Generate current time signal. + +# Parameters: +- `offset`: Offset of output signal +- `start_time`: [s] Output `y = offset` for `t < start_time` +""" +function ContinuousClock(;name, offset=0, start_time=0) + @named output = RealOutput() + pars = @parameters offset=offset start_time=start_time + eqs = [ + output.u ~ offset + ifelse(t < start_time, 0, t - start_time) + ] + compose(ODESystem(eqs, t, [], pars; name=name), [output]) +end + +""" +Generate ramp signal. + +# Parameters: +- `height`: Height of ramp +- `duration`: [s] Duration of ramp (= 0.0 gives a Step) +- `offset`: Offset of output signal +- `start_time`: [s] Output `y = offset` for `t < start_time` +""" +function Ramp(;name, + offset=0, + height=1, + duration=1, + start_time=0) + + @named output = RealOutput() + pars = @parameters offset=offset start_time=start_time height=height duration=duration + eqs = [ + output.u ~ offset + ifelse(t < start_time, 0, + ifelse(t < (start_time + duration), (t - start_time) * height / duration, height)) + ] + compose(ODESystem(eqs, t, [], pars; name=name), [output]) +end + +""" +Generate step signal. + +# Parameters: +- `height`: Height of step +- `offset`: Offset of output signal +- `start_time`: [s] Output `y = offset` for `t < start_time` +""" +function Step(;name, offset=0, height=1, start_time=0) + @named output = RealOutput() + pars = @parameters offset=offset start_time=start_time height=height + eqs = [ + output.u ~ offset + ifelse(t < start_time, 0, height) + ] + compose(ODESystem(eqs, t, [], pars; name=name), [output]) +end + +""" +Generate exponentially damped sine signal. + +# Parameters: +- `frequency`: [Hz] Frequency of sine wave +- `amplitude`: Amplitude of sine wave +- `damping`: [1/s] Damping coefficient of sine wave +- `phase`: [rad] Phase of sine wave +- `offset`: Offset of output signal +- `start_time`: [s] Output `y = offset` for `t < start_time` +""" +function ExpSine(;name, + frequency, + amplitude=1, + damping=0.1, + phase=0, + offset=0, + start_time=0) + + @named output = RealOutput() + pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase damping=damping + eqs = [ + output.u ~ offset + ifelse(t < start_time, 0, amplitude * exp(-damping * (t - start_time)) * sin(2*pi*frequency*(t - start_time) + phase)) + ] + compose(ODESystem(eqs, t, [], pars; name=name), [output]) +end + +# TODO: +# - Exponentials Generate a rising and falling exponential signal +# - Pulse Generate pulse signal of type Real +# - SawTooth Generate saw tooth signal +# - Trapezoid Generate trapezoidal signal of type Real diff --git a/src/Blocks/utils.jl b/src/Blocks/utils.jl new file mode 100644 index 000000000..97fd800d5 --- /dev/null +++ b/src/Blocks/utils.jl @@ -0,0 +1,91 @@ +@connector function RealInput(;name, nin=1, u_start=nin > 1 ? 0.0 : zeros(nin)) + if nin == 1 + @variables u(t) = u_start + else + @variables u[1:nin](t) = u_start + u = collect(u) + end + ODESystem(Equation[], t, [u...], []; name=name) +end +@doc """ + RealInput(;name, nin=1, u_start=nin > 1 ? 0.0 : zeros(nin)) + +Connector with one input signal of type Real. + +# Parameters: +- `nin`: Number of inputs +- `u_start`: Initial value for `u` + +# States: +- `u`: Value of of the connector; if nin=1 this is a scalar +""" RealInput + +@connector function RealOutput(;name, nout=1, u_start=nout > 1 ? 0.0 : zeros(nout)) + if nout == 1 + @variables u(t) = u_start + else + @variables u[1:nout](t) = u_start + u = collect(u) + end + ODESystem(Equation[], t, [u...], []; name=name) +end +@doc """ + RealOutput(;name, nout=1, u_start=nout > 1 ? 0.0 : zeros(nout)) + +Connector with one output signal of type Real. + +# Parameters: +- `nout`: Number of inputs +- `u_start`: Initial value for `u` + +# States: +- `u`: Value of of the connector; if nout=1 this is a scalar +""" RealOutput + +""" + SISO(;name, u_start=0.0, y_start=0.0) + +Single Input Single Output continuous control block. + +# Parameters: +- `u_start`: Initial value for the input +- `y_start`: Initial value for the output +""" +function SISO(;name, u_start=0.0, y_start=0.0) + @named input = RealInput(u_start=u_start) + @named output = RealOutput(u_start=y_start) + @variables begin + u(t)=u_start + y(t)=y_start + end + eqs = [ + u ~ input.u + y ~ output.u + ] + return ODESystem(eqs, t, [u, y], []; name=name, systems=[input, output]) +end + +""" + MIMO(;name, nin=1, nout=1, u_start=zeros(nin), y_start=zeros(nout)) + +Base class for a multiple Input multiple Output continuous control block. + +# Parameters: +- `nin`: Input dimension +- `nout`: Output dimension +- `u_start`: Initial value for the input +- `y_start`: Initial value for the output +""" +function MIMO(;name, nin=1, nout=1, u_start=zeros(nin), y_start=zeros(nout)) + @named input = RealInput(nin=nin, u_start=u_start) + @named output = RealOutput(nout=nout, u_start=y_start) + @variables begin + u[1:nin](t)=u_start + y[1:nout](t)=y_start + end + eqs = [ + [u[i] ~ input.u[i] for i in 1:nin]..., + [y[i] ~ output.u[i] for i in 1:nout]..., + ] + return ODESystem(eqs, t, vcat(u..., y...), []; name=name, systems=[input, output]) +end \ No newline at end of file diff --git a/src/Electrical/Analog/ideal_components.jl b/src/Electrical/Analog/ideal_components.jl index 9f5807dfd..3e8844da5 100644 --- a/src/Electrical/Analog/ideal_components.jl +++ b/src/Electrical/Analog/ideal_components.jl @@ -1,67 +1,145 @@ +""" +```julia +function Ground(;name) +``` + +Ground node with the potential of zero and connector `g`. Every circuit must have one ground +node. + +# Connectors +- `g` +""" function Ground(;name) @named g = Pin() eqs = [g.v ~ 0] - ODESystem(eqs, t, [], [], systems=[g], name=name) + ODESystem(eqs, t, [], []; systems=[g], name=name) end +""" +```julia function Resistor(;name, R = 1.0) - val = R - - @named p = Pin() - @named n = Pin() - @parameters R - @variables v(t) +``` + +Creates an ideal Resistor following Ohm's Law. + +# States +- `v(t)`: [`V`] The voltage across the resistor, given by `p.i * R` +# Connectors +- `p` + Positive pin +- `n` + Negative pin + +# Parameters: +- `R`: [`Ω`] Resistance +""" +function Resistor(;name, R=1.0) + @named oneport = OnePort() + @unpack v, i = oneport + pars = @parameters R=R eqs = [ - v ~ p.v - n.v - 0 ~ p.i + n.i - v ~ p.i * R - ] - ODESystem(eqs, t, [v], [R], systems=[p, n], defaults=Dict(R => val), name=name) + v ~ i * R + ] + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end +""" +```julia function Capacitor(; name, C = 1.0) - val = C +``` + +Creates an ideal Capacitor. + +# States +- `v(t)`: [`V`] The voltage across the capacitor, given by `D(v) ~ p.i / C` - @named p = Pin() - @named n = Pin() - @parameters C - @variables v(t) +# Connectors +- `p` + Positive pin +- `n` + Negative pin - D = Differential(t) +# Parameters: +- `C`: [`F`] Capacitance +- `v_start`: [`V`] Initial voltage of capacitor +""" +function Capacitor(;name, C=1.0, v_start=0.0) + @named oneport = OnePort(;v_start=v_start) + @unpack v, i = oneport + pars = @parameters C=C eqs = [ - v ~ p.v - n.v - 0 ~ p.i + n.i - D(v) ~ p.i / C - ] - ODESystem(eqs, t, [v], [C], systems=[p, n], defaults=Dict(C => val), name=name) + D(v) ~ i / C + ] + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end +""" +```julia function Inductor(; name, L = 1.0) - val = L +``` + +Creates an ideal Inductor. - @named p = Pin() - @named n = Pin() - @parameters L - @variables v(t) +# States +- `v(t)`: [`V`] The voltage across the inductor, given by `D(p.i) ~ v / L` - D = Differential(t) +# Connectors +- `p` + Positive pin +- `n` + Negative pin + +# Parameters: +- `L`: [`H`] Inductance +- `i_start`: [`A`] Initial current through inductor +""" +function Inductor(;name, L=1.0e-6, i_start=0.0) + @named oneport = OnePort(;i_start=i_start) + @unpack v, i = oneport + pars = @parameters L=L eqs = [ - v ~ p.v - n.v - 0 ~ p.i + n.i - D(p.i) ~ v / L - ] - ODESystem(eqs, t, [v], [L], systems=[p, n], defaults=Dict(L => val), name=name) + D(i) ~ 1 / L * v + ] + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end +""" +```julia function IdealOpAmp(; name) +``` + +Ideal operational amplifier (norator-nullator pair). +The ideal OpAmp is a two-port. The left port is fixed to `v1 = 0` and `i1 = 0` (nullator). +At the right port both any voltage `v2` and any current `i2` are possible (norator). + +# States +- `v1(t)`: [`V`] Voltage of left port +- `v2(t)`: [`V`] Voltage of right port +- `i1(t)`: [`A`] Current of left port +- `i2(t)`: [`A`] Current of right port + +# Connectors +- `p1` + Positive pin (left port) +- `p2` + Positive pin (right port) +- `n1` + Negative pin (left port) +- `n2` + Negative pin (right port) +""" +function IdealOpAmp(;name) @named p1 = Pin() @named p2 = Pin() @named n1 = Pin() @named n2 = Pin() - @variables v1(t) v2(t) # u"v" - @variables i1(t) i2(t) # u"A" - + sts = @variables begin + v1(t) + v2(t) + i1(t) + i2(t) + end eqs = [ v1 ~ p1.v - n1.v v2 ~ p2.v - n2.v @@ -72,5 +150,5 @@ function IdealOpAmp(; name) v1 ~ 0 i1 ~ 0 ] - ODESystem(eqs, t, [i1, i2, v1, v2], [], systems=[p1, p2, n1, n2], name=name) + ODESystem(eqs, t, sts, [], systems=[p1, p2, n1, n2], name=name) end diff --git a/src/Electrical/Analog/sensors.jl b/src/Electrical/Analog/sensors.jl index c69842ba5..491686eb8 100644 --- a/src/Electrical/Analog/sensors.jl +++ b/src/Electrical/Analog/sensors.jl @@ -1,37 +1,106 @@ +""" +```julia +function CurrentSensor(; name) +``` + +Creates a circuit component that measures the current flowing through it. Analogous to +an ideal ammeter. + +# States +- `i(t)`: [`A`] Current through the sensor + +# Connectors +- `p` + Positive pin +- `n` + Negative pin +""" function CurrentSensor(; name) @named p = Pin() @named n = Pin() - @variables i(t) + @variables i(t)=1.0 eqs = [ p.v ~ n.v i ~ p.i i ~ -n.i ] - ODESystem(eqs, t, [i], [], systems=[p, n], defaults=Dict(i => 1.0), name=name) + ODESystem(eqs, t, [i], [], systems=[p, n]; name=name) end +""" +```julia +function PotentialSensor(; name) +``` + +Creates a circuit component which measures the potential at a pin. + +# States +- `phi(t)`: [`V`] The potential at this point + +# Connectors +- `p` + Pin at which potential is to be measured +""" function PotentialSensor(; name) @named p = Pin() - @variables phi(t) + @variables phi(t)=1.0 eqs = [ p.i ~ 0 phi ~ p.v ] - ODESystem(eqs, t, [phi], [], systems=[p], defaults=Dict(phi => 1.0), name=name) + ODESystem(eqs, t, [phi], [], systems=[p]; name=name) end +""" +```julia +function VoltageSensor(; name) +``` + +Creates a circuit component that measures the voltage across it. Analogous to +an ideal voltmeter. + +# States +- `v(t)`: [`V`] The voltage across this component + +# Connectors +- `p` + Positive pin +- `n` + Negative pin +""" function VoltageSensor(; name) @named p = Pin() @named n = Pin() - @variables v(t) + @variables v(t)=1.0 eqs = [ p.i ~ 0 n.i ~ 0 v ~ p.v - n.v ] - ODESystem(eqs, t, [v], [], systems=[p, n], defaults=Dict(v => 1.0), name=name) + ODESystem(eqs, t, [v], []; systems=[p, n], name=name) end +""" +```julia +function PowerSensor(; name) +``` + +Combines a [`VoltageSensor`](@ref) and a [`CurrentSensor`](@ref) to measure the power being +consumed by a circuit. + +# States +- `power(t)`: [`W`] The power being consumed, given by the product of voltage and current. + +# Connectors +- `pc` + Corresponds to the `p` pin of the [`CurrentSensor`](@ref) +- `nc` + Corresponds to the `n` pin of the [`CurrentSensor`](@ref) +- `pv` + Corresponds to the `p` pin of the [`VoltageSensor`](@ref) +- `nv` + Corresponds to the `n` pin of the [`VoltageSensor`](@ref) +""" function PowerSensor(; name) @named pc = Pin() @named nc = Pin() @@ -39,7 +108,7 @@ function PowerSensor(; name) @named nv = Pin() @named voltage_sensor = VoltageSensor() @named current_sensor = CurrentSensor() - @variables power(t) + @variables power(t)=1.0 eqs = [ connect(voltage_sensor.p, pv) connect(voltage_sensor.n, nv) @@ -47,9 +116,30 @@ function PowerSensor(; name) connect(current_sensor.n, nc) power ~ current_sensor.i * voltage_sensor.v ] - ODESystem(eqs, t, [power], [], systems=[pc, nc, pv, nv, voltage_sensor, current_sensor], defaults=Dict(power => 1.0), name=name) + ODESystem(eqs, t, [power], []; systems=[pc, nc, pv, nv, voltage_sensor, current_sensor], name=name) end +""" +```julia +function MultiSensor(; name) +``` + +Combines a [`VoltageSensor`](@ref) and a [`CurrentSensor`](@ref). + +# States +- `v(t)`: [`V`] The voltage across the [`VoltageSensor`](@ref) +- `i(t)`: [`V`] The current across the [`CurrentSensor`](@ref) + +# Connectors +- `pc` + Corresponds to the `p` pin of the [`CurrentSensor`](@ref) +- `nc` + Corresponds to the `n` pin of the [`CurrentSensor`](@ref) +- `pv` + Corresponds to the `p` pin of the [`VoltageSensor`](@ref) +- `nv` + Corresponds to the `n` pin of the [`VoltageSensor`](@ref) +""" function MultiSensor(; name) @named pc = Pin() @named nc = Pin() @@ -57,7 +147,10 @@ function MultiSensor(; name) @named nv = Pin() @named voltage_sensor = VoltageSensor() @named current_sensor = CurrentSensor() - @variables i(t) v(t) + sts = @variables begin + i(t)=1.0 + v(t)=1.0 + end eqs = [ connect(voltage_sensor.p, pv) connect(voltage_sensor.n, nv) @@ -66,5 +159,5 @@ function MultiSensor(; name) i ~ current_sensor.i v ~ voltage_sensor.v ] - ODESystem(eqs, t, [i, v], [], systems=[pc, nc, pv, nv, voltage_sensor, current_sensor], defaults=Dict(i => 1.0, v => 1.0), name=name) + ODESystem(eqs, t, sts, []; systems=[pc, nc, pv, nv, voltage_sensor, current_sensor], name=name) end diff --git a/src/Electrical/Analog/sources.jl b/src/Electrical/Analog/sources.jl index 4d1abbc16..5305ff97c 100644 --- a/src/Electrical/Analog/sources.jl +++ b/src/Electrical/Analog/sources.jl @@ -1,299 +1,584 @@ # Define and register smooth functions -_cos_wave(x, f, A, st, ϕ) = A*cos(2*π*f*(x-st) + ϕ) -_damped_sine_wave(x, f, A, st, ϕ, d) = exp((st-x)*d)*A*sin(2*π*f*(x-st) + ϕ) -_ramp(x, δ, st, et, h) = h/(et-st)*(_xH(x, δ, st) - _xH(x, δ, et)) -_square_wave(x, δ, f, A, st) = A*2atan(sin(2π*(x-st)*f)/δ)/π -_step(x, δ, h, a) = h*(atan((x-a)/δ)/π + 0.5) -_triangular_wave(x, δ, f, A, st) = A*(1-2acos((1 - δ)sin(2π*(x-st)*f))/π) -_xH(x, δ, tₒ) = 0.5*(x-tₒ)*(1+((x-tₒ)/sqrt((x-tₒ)^2+δ^2))) - -@register _cos_wave(x, f, A, st, ϕ) -@register _damped_sine_wave(x, f, A, st, ϕ, damping) -@register _ramp(x, δ, st, et, h) -@register _square_wave(x, δ, f, A, st) -@register _step(x, δ, h, a) -@register _triangular_wave(x, δ, f, A, st) +_cos_wave(t, f, A, st, ϕ) = A*cos(2*π*f*(t - st) + ϕ) +_sin_wave(t, f, A, st, ϕ) = A*sin(2*π*f*(t - st) + ϕ) +_damped_sine_wave(t, f, A, st, ϕ, d) = exp((st-t)*d)*A*sin(2*π*f*(t-st) + ϕ) +_ramp(t, δ, st, et, h) = h/(et-st)*(_xH(t, δ, st) - _xH(t, δ, et)) +_square_wave(t, δ, f, A, st) = A*2atan(sin(2π*(t-st)*f)/δ)/π +_step(t, δ, h, a) = h*(atan((t-a)/δ)/π + 0.5) +_triangular_wave(t, δ, f, A, st) = A*(1-2acos((1 - δ)sin(2π*(t-st)*f))/π) +_xH(t, δ, tₒ) = (t-tₒ)*(1+((t-tₒ)/sqrt((t-tₒ)^2+δ^2)))/2 + +@register_symbolic _cos_wave(t, f, A, st, ϕ) +@register_symbolic _sin_wave(t, f, A, st, ϕ) +@register_symbolic _damped_sine_wave(t, f, A, st, ϕ, damping) +@register_symbolic _ramp(t, δ, st, et, h) +@register_symbolic _square_wave(t, δ, f, A, st) +@register_symbolic _step(t, δ, h, a) +@register_symbolic _triangular_wave(t, δ, f, A, st) # Voltage sources +""" +```julia function ConstantVoltage(;name, V=1.0) - val = V - - @named p = Pin() - @named n = Pin() - @parameters V - @variables v(t) - +``` + +The source for an ideal constant voltage. + +# States +- `v(t)`: [`V`] The voltage across this source, given by `p.v - n.v` and is always constant + +# Connectors +- `p` + Positive pin +- `n` + Negative pin + +# Parameters: +- `V`: [`V`] The constant voltage across the terminals of this source +""" +function ConstantVoltage(;name, V = 1.0) + @named oneport = OnePort() + @unpack v, i = oneport + pars = @parameters V=V eqs = [ - v ~ p.v - n.v - 0 ~ p.i + n.i - v ~ V - ] - ODESystem(eqs, t, [v], [V], systems=[p, n], defaults=Dict(V => val), name=name) + v ~ V + ] + + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end +""" +```julia function CosineVoltage(;name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0, phase=0.0) - o, A, f, st, ϕ = offset, amplitude, frequency, starttime, phase +``` + +A source in which the voltage across its terminals is a cosine function of time. + + +# States +- `v(t)` + The voltage across this source, given by `p.v - n.v` + +# Connectors +- `p` + Positive port +- `n` + Negative port + +# Observables +- `offset`: [`V`] + A constant offset added to the voltage output +- `amplitude`: [`V`] + The amplitude of the cosine function +- `frequency`: [`Hz`] + The frequency of the cosine function +- `starttime`: [`s`] + The time at which the source starts functioning. Before this time, the voltage across + its terminals is 0. +- `phase`: [`rad`] + The phase offset of the cosine function +""" +function CosineVoltage(;name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0) δ = 0.00001 - @named p = Pin() - @named n = Pin() - @parameters offset amplitude frequency starttime phase - @variables v(t) - + @named oneport = OnePort() + @unpack v, i = oneport + pars = @parameters begin + offset=offset + amplitude=amplitude + frequency=frequency + start_time=start_time + phase=phase + end eqs = [ - v ~ p.v - n.v - v ~ _cos_wave(t, f, A, st, ϕ) * _step(t, δ, 1.0, st) + offset - 0 ~ p.i + n.i - ] - defaults = Dict(zip((offset, amplitude, frequency, starttime, phase), (o, A, f, st, ϕ))) - ODESystem(eqs, t, [v], [offset, amplitude, frequency, starttime, phase], systems=[p, n], defaults=defaults, name=name) + v ~ _cos_wave(t, frequency, amplitude, start_time, phase) * _step(t, δ, 1.0, start_time) + offset + ] + + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end -function DampedSineVoltage(;name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0, phase=0.0, damping_coef=0.0) - o, A, f, st, ϕ, d = offset, amplitude, frequency, starttime, phase, damping_coef - δ = 0.0001 - - @named p = Pin() - @named n = Pin() - @parameters offset amplitude frequency starttime phase damping_coef - @variables v(t) - +""" +```julia +function ExpSineVoltage(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0, damping=0.0) +``` + +A source in which the voltage across its terminals is a damped sine function of time. + +# States +- `v(t)`: [`V`] + The voltage across this source, given by `p.v - n.v` + +# Connectors +- `p` + Positive port +- `n` + Negative port + +# Parameters +- `offset`: [`V`] + A constant offset added to the voltage output +- `amplitude`: [`V`] + The amplitude of the damped sine function +- `frequency`: [`Hz`] + The frequency of the damped sine function +- `start_time`: [`s`] + The time at which the source starts functioning. Before this time, the voltage across + its terminals is `offset`. +- `phase`: [`rad`] + The phase offset of the damped sine function +- `damping_coef`: [`1/s`] + Damping coefficient of the damped sine function +""" +function ExpSineVoltage(;name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0, damping=0.0) + δ = 0.00001 + @named oneport = OnePort() + @unpack v, i = oneport + pars = @parameters begin + offset=offset + amplitude=amplitude + frequency=frequency + start_time=start_time + phase=phase + damping=damping + end eqs = [ - v ~ p.v - n.v - v ~ _step(t, δ, o, 0.0) + _damped_sine_wave(t, f, A, st, ϕ, d) * _step(t, δ, 1.0, st) - 0 ~ p.i + n.i - ] - defaults = Dict(zip((offset, amplitude, frequency, starttime, phase, damping_coef), (o, A, f, st, ϕ, d))) - ODESystem(eqs, t, [v], [offset, amplitude, frequency, starttime, phase, damping_coef], systems=[p, n], defaults=defaults, name=name) + v ~ _damped_sine_wave(t, frequency, amplitude, start_time, phase, damping) * _step(t, δ, 1.0, start_time) + ] + + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end -function RampVoltage(;name, offset=0.0, starttime=0.0, endtime=1.0, height=1.0) - o, st, et, h = offset, starttime, endtime, height - δ = 0.0001 - - @named p = Pin() - @named n = Pin() - @parameters offset starttime endtime height - @variables v(t) - +""" +```julia +function RampVoltage(;name, offset=0.0, start_time=0.0, duration=1.0, height=1.0) +``` + +A source in which the voltage across grows linearly from `offset` to `offset+height` over +the time interval `duration` starting at `start_time` + +# States +- `v(t)`: [`V`] + The voltage across this source, given by `p.v - n.v` + +# Connectors +- `p` + Positive port +- `n` + Negative port + RampVoltage(;name, offset=0.0, start_time=0.0, duration=1.0, height=1.0) + +# Parameters +- `offset`: [`V`] + A constant offset added to the voltage output +- `start_time`: [`s`] + The time at which the voltage starts growing +- `duration`: [`s`] + The duration of the ramp (`0.0` gives a step) +- `height`: [`V`] + The amount that the voltage grows in the time interval +""" +function RampVoltage(;name, offset=0.0, start_time=0.0, duration=1.0, height=1.0) + δ = 0.00001 + @named oneport = OnePort() + @unpack v, i = oneport + pars = @parameters begin + offset=offset + height=height + start_time=start_time + duration=duration + end eqs = [ - v ~ p.v - n.v - v ~ offset + _ramp(t, δ, st, et, h) - 0 ~ p.i + n.i - ] - defaults = Dict(zip((offset, starttime, endtime, height), (o, st, et, h))) - ODESystem(eqs, t, [v], [offset, starttime, endtime, height], systems=[p, n], defaults=defaults, name=name) + v ~ _ramp(t, δ, start_time, start_time + duration, height) + offset + ] + + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end -function SineVoltage(;name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0, phase=0.0) - o, A, f, st, ϕ = offset, amplitude, frequency, starttime, phase - - @named p = Pin() - @named n = Pin() - @parameters offset amplitude frequency starttime phase - @variables v(t) +""" +```julia +function SineVoltage(;name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0) +``` + +A source in which the voltage across its terminals is a sine function of time. + + + +# States +- `v(t)`: [`V`] + The voltage across this source, given by `p.v - n.v` + +# Connectors +- `p` + Positive port +- `n` + Negative port + +# Parameters +- `offset`: [`V`] + A constant offset added to the voltage output +- `amplitude`: [`V`] + The amplitude of the sine function +- `frequency`: [`Hz`] + The frequency of the sine function +- `start_time`: [`s`] + The time at which the source starts functioning. Before this time, the voltage across + its terminals is `offset`. +- `phase`: [`rad`] + The phase offset of the sine function +""" +function SineVoltage(;name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0) + δ = 0.00001 + @named oneport = OnePort() + @unpack v, i = oneport + pars = @parameters begin + offset=offset + amplitude=amplitude + frequency=frequency + start_time=start_time + phase=phase + end eqs = [ - v ~ p.v - n.v - v ~ offset + (t > st) * (A*sin(2*π*f*(t - st) + ϕ)) - 0 ~ p.i + n.i - ] - defaults = Dict(zip((offset, amplitude, frequency, starttime, phase), (o, A, f, st, ϕ))) - ODESystem(eqs, t, [v], [offset, amplitude, frequency, starttime, phase], systems=[p, n], defaults=defaults, name=name) + v ~ _sin_wave(t, frequency, amplitude, start_time, phase) * _step(t, δ, 1.0, start_time) + offset + ] + + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end -function SquareVoltage(; name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0) - o, A, f, st = offset, amplitude, frequency, starttime +""" +```julia +function SquareVoltage(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0) +``` + +A source in which the voltage across its terminals is a square function of time. + +# States +- `v(t)`: [`V`] + The voltage across this source, given by `p.v - n.v` + +# Connectors +- `p` + Positive port +- `n` + Negative port + +# Parameters +- `offset`: [`V`] + A constant offset added to the voltage output +- `amplitude`: [`V`] + The amplitude of the square wave function +- `frequency`: [`Hz`] + The frequency of the square wave function +- `start_time`: [`s`] + The time at which the source starts functioning. Before this time, the voltage across + its terminals is `offset`. +""" +function SquareVoltage(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0) δ = 0.0001 - @named p = Pin() - @named n = Pin() - @parameters offset amplitude frequency starttime - @variables v(t) - + @named oneport = OnePort() + @unpack v, i = oneport + pars = @parameters begin + offset=offset + amplitude=amplitude + start_time=start_time + end eqs = [ - v ~ p.v - n.v - 0 ~ p.i + n.i - v ~ o + _square_wave(t, δ, f, A, st) * (t > st) - ] - defaults = Dict(zip((offset, amplitude, frequency, starttime), (o, A, f, st))) - ODESystem(eqs, t, [v], [offset, amplitude, frequency, starttime], systems=[p, n], defaults=defaults, name=name) + v ~ _square_wave(t, δ, frequency, amplitude, start_time) * _step(t, δ, 1.0, start_time) + offset + ] + + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end -function StepVoltage(;name, offset=0.0, starttime=0.0, height=1.0) - o, st, h = offset, starttime, height +""" +```julia +function StepVoltage(;name, offset=0.0, start_time=0.0, height=1.0) +``` + +A source in which the voltage across its terminals increases from `offset` to `offset+height` at +`starttime` + +# States +- `v(t)`: [`V`] + The voltage across this source, given by `p.v - n.v` + +# Connectors +- `p` + Positive port +- `n` + Negative port + +# Observables +- `offset`: [`V`] + A constant offset added to the voltage output +- `start_time`: [`s`] + The time at which the source starts functioning, and the voltage jumps +- `height`: [`V`] + Magnitude of increase in voltage +""" +function StepVoltage(;name, offset=0.0, start_time=0.0, height=1.0) δ = 0.0001 - @named p = Pin() - @named n = Pin() - @parameters offset starttime height - @variables v(t) - + @named oneport = OnePort() + @unpack v, i = oneport + pars = @parameters begin + offset=offset + height=height + start_time=start_time + end eqs = [ - v ~ p.v - n.v - v ~ offset + _step(t, δ, h, st) - 0 ~ p.i + n.i - ] - defaults = Dict(zip((offset, starttime, height), (o, st, h))) - ODESystem(eqs, t, [v], [offset, starttime, height], systems=[p, n], defaults=defaults, name=name) + v ~ _step(t, δ, height, start_time) + offset + ] + + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end -function TriangularVoltage(; name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0) - o, A, f, st = offset, amplitude, frequency, starttime - δ = 0.0001 +""" +```julia +function TriangularVoltage(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0) +``` + +A source in which the voltage across its terminals is a triangular function of time. + +# States +- `v(t)`: [`V`] + The voltage across this source, given by `p.v - n.v` + +# Connectors +- `p` + Positive port +- `n` + Negative port + +# Observables +- `offset`: [`V`] + A constant offset added to the voltage output +- `amplitude`: [`V`] + Amplitude of the triangular wave function +- `frequency`: [`Hz`] + Frequency of the triangular wave function +- `start_time`: [`s`] + The time at which the source starts functioning. Before this, the output of the source is + `offset` +""" +function TriangularVoltage(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0) + δ = 0.00001 - @named p = Pin() - @named n = Pin() - @parameters offset amplitude frequency starttime - @variables v(t) - + @named oneport = OnePort() + @unpack v, i = oneport + pars = @parameters begin + offset=offset + amplitude=amplitude + frequency=frequency + start_time=start_time + end eqs = [ - v ~ p.v - n.v - 0 ~ p.i + n.i - v ~ offset + (t>st) * _triangular_wave(t, δ, f, A, st) + v ~ _triangular_wave(t, δ, frequency, amplitude, start_time) * _step(t, δ, 1.0, start_time) + offset ] - defaults = Dict(zip((offset, amplitude, frequency, starttime), (o, A, f, st))) - ODESystem(eqs, t, [v], [offset, amplitude, frequency, starttime], systems=[p, n], defaults=defaults, name=name) + + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end -# Current Sources -function ConstantCurrent(;name, I=1.0) - val = I +# Current Sources ###################################################################################################### +""" + ConstantCurrent(;name, I = 1.0) - @named p = Pin() - @named n = Pin() - @parameters I - @variables i(t) +Source for constant current. +# Parameters: +- `I`: [A] Current +""" +function ConstantCurrent(;name, I = 1.0) + @named oneport = OnePort() + @unpack v, i = oneport + pars = @parameters I=I eqs = [ - 0 ~ p.i + n.i - i ~ p.i - i ~ I - ] - ODESystem(eqs, t, [i], [I], systems=[p, n], defaults=Dict(I => val), name=name) + i ~ I + ] + + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end -function CosineCurrent(;name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0, phase=0.0) - o, A, f, st, ϕ = offset, amplitude, frequency, starttime, phase - δ = 0.00001 +""" + CosineCurrent(;name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0) - @named p = Pin() - @named n = Pin() - @parameters offset amplitude frequency starttime phase - @variables i(t) +Generate cosine current. +# Parameters: +- `frequency`: [Hz] Frequency of sine wave +- `amplitude`: [A] Amplitude of sine wave +- `phase`: [rad] Phase of sine wave +- `offset`: [A] Offset of output current +- `start_time`: [s] Output `y = offset` for `t < start_time` +""" +function CosineCurrent(;name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0) + δ = 0.00001 + + @named oneport = OnePort() + @unpack v, i = oneport + pars = @parameters begin + offset=offset + amplitude=amplitude + frequency=frequency + start_time=start_time + phase=phase + end eqs = [ - i ~ _cos_wave(t, f, A, st, ϕ) * _step(t, δ, 1.0, st) + offset - 0 ~ p.i + n.i - i ~ p.i - ] - defaults = Dict(zip((offset, amplitude, frequency, starttime, phase), (o, A, f, st, ϕ))) - ODESystem(eqs, t, [i], [offset, amplitude, frequency, starttime, phase], systems=[p, n], defaults=defaults, name=name) + i ~ _cos_wave(t, frequency, amplitude, start_time, phase) * _step(t, δ, 1.0, start_time) + offset + ] + + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end -function DampedSineCurrent(;name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0, phase=0.0, damping_coef=1.0) - o, A, f, st, ϕ, d = offset, amplitude, frequency, starttime, phase, damping_coef - δ = 0.0001 +""" + ExpSineCurrent(;name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0, damping=0.0) - @named p = Pin() - @named n = Pin() - @parameters offset amplitude frequency starttime phase damping_coef - @variables i(t) +Generate damped sine current. +# Parameters: +- `frequency`: [Hz] Frequency of sine wave +- `amplitude`: [A] Amplitude of sine wave +- `damping`: [1/s] Damping coefficient of sine wave +- `phase`: [rad] Phase of sine wave +- `offset`: [A] Offset of output current +- `start_time`: [s] Output `y = offset` for `t < start_time` +""" +function ExpSineCurrent(;name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0, damping=0.0) + δ = 0.00001 + @named oneport = OnePort() + @unpack v, i = oneport + pars = @parameters begin + offset=offset + amplitude=amplitude + frequency=frequency + start_time=start_time + phase=phase + damping=damping + end eqs = [ - i ~ _step(t, δ, o, 0.0) + _damped_sine_wave(t, f, A, st, ϕ, d) * _step(t, δ, 1.0, st) - 0 ~ p.i + n.i - i ~ p.i - ] - defaults = Dict(zip((offset, amplitude, frequency, starttime, phase, damping_coef), (o, A, f, st, ϕ, d))) - ODESystem(eqs, t, [i], [offset, amplitude, frequency, starttime, phase, damping_coef], systems=[p, n], defaults=defaults, name=name) + i ~ _damped_sine_wave(t, frequency, amplitude, start_time, phase, damping) * _step(t, δ, 1.0, start_time) + ] + + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end -function RampCurrent(;name, offset=0.0, starttime=0.0, endtime=1.0, height=1.0) - o, st, et, h = offset, starttime, endtime, height - δ = 0.0001 +""" + RampCurrent(;name, offset=0.0, start_time=0.0, duration=1.0, height=1.0) - @named p = Pin() - @named n = Pin() - @parameters offset starttime endtime height - @variables i(t) +Generate ramp current. +# Parameters: +- `height`: [A] Height of ramp +- `duration`: [s] Duration of ramp (= 0.0 gives a Step) +- `offset`: [A] Offset of output current +- `start_time`: [s] Output `y = offset` for `t < start_time` +""" +function RampCurrent(;name, offset=0.0, start_time=0.0, duration=1.0, height=1.0) + δ = 0.00001 + @named oneport = OnePort() + @unpack v, i = oneport + pars = @parameters begin + offset=offset + height=height + start_time=start_time + duration=duration + end eqs = [ - i ~ _step(t, δ, o, 0.0) + _ramp(t, δ, st, et, h) - 0 ~ p.i + n.i - i ~ p.i - ] - defaults = Dict(zip((offset, starttime, endtime, height), (o, st, et, h))) - ODESystem(eqs, t, [i], [offset, starttime, endtime, height], systems=[p, n], defaults=defaults, name=name) + i ~ _ramp(t, δ, start_time, start_time + duration, height) + offset + ] + + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end -function SineCurrent(;name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0, phase=0.0) - o, A, f, st, ϕ = offset, amplitude, frequency, starttime, phase +""" + SineCurrent(;name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0) - @named p = Pin() - @named n = Pin() - @parameters offset amplitude frequency starttime phase - @variables i(t) +Generate sine current. +# Parameters: +- `frequency`: [Hz] Frequency of sine wave +- `amplitude`: [A] Amplitude of sine wave +- `phase`: [rad] Phase of sine wave +- `offset`: [A] Offset of output current +- `start_time`: [s] Output `y = offset` for `t < start_time` +""" +function SineCurrent(;name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0, phase=0.0) + δ = 0.00001 + + @named oneport = OnePort() + @unpack v, i = oneport + pars = @parameters begin + offset=offset + amplitude=amplitude + frequency=frequency + start_time=start_time + phase=phase + end eqs = [ - i ~ offset + (t > st) * (A*sin(2*π*f*(t - st) + ϕ)) - 0 ~ p.i + n.i - i ~ p.i - ] - defaults = Dict(zip((offset, amplitude, frequency, starttime, phase), (o, A, f, st, ϕ))) - ODESystem(eqs, t, [i], [offset, amplitude, frequency, starttime, phase], systems=[p, n], defaults=defaults, name=name) + i ~ _sin_wave(t, frequency, amplitude, start_time, phase) * _step(t, δ, 1.0, start_time) + offset + ] + + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end -function SquareCurrent(; name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0) - o, A, f, st = offset, amplitude, frequency, starttime +function SquareCurrent(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0) δ = 0.0001 - @named p = Pin() - @named n = Pin() - @parameters offset amplitude frequency starttime - @variables i(t) - + @named oneport = OnePort() + @unpack v, i = oneport + pars = @parameters begin + offset=offset + amplitude=amplitude + start_time=start_time + end eqs = [ - 0 ~ p.i + n.i - i ~ o + _square_wave(t, δ, f, A, st) * (t > st) - i ~ p.i - ] - defaults = Dict(zip((offset, amplitude, frequency, starttime), (o, A, f, st))) - ODESystem(eqs, t, [i], [offset, amplitude, frequency, starttime], systems=[p, n], defaults=defaults, name=name) + i ~ _square_wave(t, δ, frequency, amplitude, start_time) * _step(t, δ, 1.0, start_time) + offset + ] + + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end -function StepCurrent(;name, offset=0.0, starttime=0.0, height=1.0) - o, st, h = offset, starttime, height - δ = 0.0001 +""" + StepCurrent(;name, offset=0.0, start_time=0.0, height=1.0) - @named p = Pin() - @named n = Pin() - @parameters offset starttime height - @variables i(t) +Generate step current. +# Parameters: +- `height`: [A] Height of step +- `offset`: [A] Offset of output current +- `start_time`: [s] Output `y = offset` for `t < start_time` +""" +function StepCurrent(;name, offset=0.0, start_time=0.0, height=1.0) + δ = 0.0001 + + @named oneport = OnePort() + @unpack v, i = oneport + pars = @parameters begin + offset=offset + height=height + start_time=start_time + end eqs = [ - i ~ offset + _step(t, δ, h, st) - 0 ~ p.i + n.i - i ~ p.i - ] - defaults = Dict(zip((offset, starttime, height), (o, st, h))) - ODESystem(eqs, t, [i], [offset, starttime, height], systems=[p, n], defaults=defaults, name=name) + i ~ _step(t, δ, height, start_time) + offset + ] + + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end -function TriangularCurrent(; name, offset=0.0, amplitude=1.0, frequency=1.0, starttime=0.0) - o, A, f, st = offset, amplitude, frequency, starttime - δ = 0.0001 +function TriangularCurrent(; name, offset=0.0, amplitude=1.0, frequency=1.0, start_time=0.0) + δ = 0.00001 - @named p = Pin() - @named n = Pin() - @parameters offset amplitude frequency starttime - @variables i(t) - + @named oneport = OnePort() + @unpack v, i = oneport + pars = @parameters begin + offset=offset + amplitude=amplitude + frequency=frequency + start_time=start_time + end eqs = [ - 0 ~ p.i + n.i - i ~ offset + _step(t, δ, 1, st) * _triangular_wave(t, δ, f, A, st) - i ~ p.i + i ~ _triangular_wave(t, δ, frequency, amplitude, start_time) * _step(t, δ, 1.0, start_time) + offset ] - defaults = Dict(zip((offset, amplitude, frequency, starttime), (o, A, f, st))) - ODESystem(eqs, t, [i], [offset, amplitude, frequency, starttime], systems=[p, n], defaults=defaults, name=name) + + extend(ODESystem(eqs, t, [], pars; name=name), oneport) end diff --git a/src/Electrical/Digital/components.jl b/src/Electrical/Digital/components.jl index 1e40ae607..6ae1b1170 100644 --- a/src/Electrical/Digital/components.jl +++ b/src/Electrical/Digital/components.jl @@ -1,5 +1,25 @@ # Adders +""" +```julia +function HalfAdder(; name) +``` + +Takes two bits as input, and outputs the sum and the carry + +# States +- `sum(t)` + The sum of the input bits +- `carry(t)` + The carry generated by the input bits +# Connectors +- `x1`, `x2` + The two inputs to add +- `y1` + Output [`DigitalPin`](@ref) corresponding to the sum +- `y2` + Output [`DigitalPin`](@ref) corresponding to the carry +""" function HalfAdder(; name) @named x1 = DigitalPin() @named x2 = DigitalPin() @@ -16,6 +36,27 @@ function HalfAdder(; name) ODESystem(eqs, t, [sum, carry], [], systems=[x1, x2, y1, y2], name=name) end +""" +```julia +function FullAdder(; name) +``` + +Takes three bits as input, and outputs the sum and the carry + +# States +- `sum(t)` + The sum of the input bits +- `carry(t)` + The carry generated by the input bits + +# Connectors +- `x1`, `x2`, `x3` + The three inputs to add +- `y1` + Output [`DigitalPin`](@ref) corresponding to the sum +- `y2` + Output [`DigitalPin`](@ref) corresponding to the carry +""" function FullAdder(; name) @named x1 = DigitalPin() @named x2 = DigitalPin() @@ -37,6 +78,24 @@ end # This selects data from the `N` input ports (`d₀` to `dₙ₋₁`) # using values of `n` select lines, where `N = 2^n` +""" +```julia +function MUX(; name, N=4) +``` + +Standard Multiplexer. Selects data from `N` input ports using the values +of `n` select lines, where `N=2ⁿ`. For the `i`th input port to be selected, +the values of the select lines should correspond to the binary representation +of `i`. + +# Connectors +- `d1`, `d2`, ... + The `N` input lines +- `s1`, `s2`, ... + The `n` select lines +- `y` + The output, selected from one of the `N` input lines +""" function MUX(; name, N=4) n = log2(N) try n = Int(n) catch(e) throw("`N` must be a power of 2") end @@ -68,6 +127,25 @@ end # This selects one of the `N` output ports (`y₀` to `yₙ₋₁`) # to transmit data `d` using values of `n` select lines, where `N = 2^n` +""" +```julia +function DEMUX(; name, N=4) +``` + +Standard Demultiplexer. Performs the reverse operation of a [`MUX`](@ref). +Selects one of the `N` output ports to transmit the input `d` using the +values of `n` select lines, where `N=2ⁿ`. For the `i`th output port to be +selected, the values of the select lines should correspond to the binary +representation of `i`. + +# Connectors +- `d` + The input to be transmitted to one of the output lines +- `s1`, `s2`, ... + The `n` select lines +- `y1`, `y2`, ... + The `N` output lines +""" function DEMUX(; name, N=4) n = log2(N) try n = Int(n) catch(e) throw("`N` must be a power of 2") end @@ -96,6 +174,20 @@ end # Encoder-Decoder # Encodes `N` inputs to `n` outputs, where `N = 2^n` +""" +```julia +function Encoder(; name, N=4) +``` + +Encodes `N` inputs to `n` outputs, where `N=2ⁿ`. Exactly one of the inputs should be `1`. +If the `i`th input is `1`, then the output corresponds to the binary representation of `i`. + +# Connectors +- `d1`, `d2`, ... + The `N` input lines +- `y1`, `y2`, ... + The `n` output lines +""" function Encoder(; name, N=4) n = log2(N) try n = Int(n) catch(e) throw("`N` must be a power of 2") end @@ -132,6 +224,21 @@ function Encoder(; name, N=4) end # Decodes `n` inputs to `N` outputs, where `N = 2^n` +""" +```julia +function Decoder(; name, n=2) +``` + +Performs the reverse operation of an [`Encoder`](@ref). Decodes `n` inputs +to `N` outputs, where `N=2ⁿ`. The `i`th output is `1` if the values of +the select lines correspond to the binary representation of `1`. + +# Connectors +- `d1`, `d2`, ... + The `n` input lines +- `y1`, `y2`, ... + The `N` output lines +""" function Decoder(; name, n=2) N = 2^n d = map(0:n-1) do i diff --git a/src/Electrical/Digital/gates.jl b/src/Electrical/Digital/gates.jl index 9d69d6291..ef480138a 100644 --- a/src/Electrical/Digital/gates.jl +++ b/src/Electrical/Digital/gates.jl @@ -1,3 +1,16 @@ +""" +```julia +function Not(; name) +``` + +NOT gate in 9-level logic. + +# Connectors +- `x` + Input [`DigitalPin`](@ref) +- `y` + Output [`DigitalPin`](@ref) +""" function Not(; name) @named x = DigitalPin() @named y = DigitalPin() @@ -9,6 +22,19 @@ function Not(; name) ODESystem(eqs, t, [], [], systems=[x, y], name=name) end +""" +```julia +function And(; name, N=2) +``` + +AND gate in 9-level logic, with `N` inputs + +# Connectors +- `x1`, `x2`, ... + `N` input [`DigitalPin`](@ref)s +- `y` + Output [`DigitalPin`](@ref) +""" function And(; name, N=2) x = map(1:N) do i DigitalPin(name=Symbol(:x, i)) @@ -23,6 +49,19 @@ function And(; name, N=2) ODESystem(eqs, t, [], [], systems=[x..., y], name=name) end +""" +```julia +function Nand(; name, N=2) +``` + +NAND gate in 9-level logic, with `N` inputs + +# Connectors +- `x1`, `x2`, ... + `N` input [`DigitalPin`](@ref)s +- `y` + Output [`DigitalPin`](@ref) +""" function Nand(; name, N=2) x = map(1:N) do i DigitalPin(name=Symbol(:x, i)) @@ -37,6 +76,19 @@ function Nand(; name, N=2) ODESystem(eqs, t, [], [], systems=[x..., y], name=name) end +""" +```julia +function Or(; name, N=2) +``` + +OR gate in 9-level logic, with `N` inputs + +# Connectors +- `x1`, `x2`, ... + `N` input [`DigitalPin`](@ref)s +- `y` + Output [`DigitalPin`](@ref) +""" function Or(; name, N=2) x = map(1:N) do i DigitalPin(name=Symbol(:x, i)) @@ -51,6 +103,19 @@ function Or(; name, N=2) ODESystem(eqs, t, [], [], systems=[x..., y], name=name) end +""" +```julia +function Nor(; name, N=2) +``` + +NOR gate in 9-level logic, with `N` inputs + +# Connectors +- `x1`, `x2`, ... + `N` input [`DigitalPin`](@ref)s +- `y` + Output [`DigitalPin`](@ref) +""" function Nor(; name, N=2) x = map(1:N) do i DigitalPin(name=Symbol(:x, i)) @@ -65,6 +130,19 @@ function Nor(; name, N=2) ODESystem(eqs, t, [], [], systems=[x..., y], name=name) end +""" +```julia +function Xor(; name, N=2) +``` + +XOR gate in 9-level logic, with `N` inputs + +# Connectors +- `x1`, `x2`, ... + `N` input [`DigitalPin`](@ref)s +- `y` + Output [`DigitalPin`](@ref) +""" function Xor(; name, N=2) x = map(1:N) do i DigitalPin(name=Symbol(:x, i)) @@ -79,6 +157,19 @@ function Xor(; name, N=2) ODESystem(eqs, t, [], [], systems=[x..., y], name=name) end +""" +```julia +function Xnor(; name, N=2) +``` + +XNOR gate in 9-level logic, with `N` inputs + +# Connectors +- `x1`, `x2`, ... + `N` input [`DigitalPin`](@ref)s +- `y` + Output [`DigitalPin`](@ref) +""" function Xnor(; name, N=2) x = map(1:N) do i DigitalPin(name=Symbol(:x, i)) diff --git a/src/Electrical/Digital/sources.jl b/src/Electrical/Digital/sources.jl index 5b25c7e2c..fe5383cc1 100644 --- a/src/Electrical/Digital/sources.jl +++ b/src/Electrical/Digital/sources.jl @@ -1,3 +1,16 @@ +""" +```julia +function PulseDiff(; name, Val=1, dt=0.1) +``` + +# States +- `val(t)` + Output value of the source + +# Connectors +- `d` + Output [`DigitalPin`](@ref) +""" function PulseDiff(; name, Val=1, dt=0.1) @named d = DigitalPin() @variables val(t) @@ -11,6 +24,17 @@ function PulseDiff(; name, Val=1, dt=0.1) ODESystem(eqs, t, [val], [], systems=[d], defaults=Dict(Val=>0), name=name) end +""" +```julia +function Set(; name) +``` + +Source that outputs a constant signal of `1`. + +# Connectors +- `d` + Output [`DigitalPin`](@ref) +""" function Set(; name) @named d = DigitalPin() @@ -20,6 +44,17 @@ function Set(; name) ODESystem(eqs, t, [], [], systems=[d], name=name) end +""" +```julia +function Reset(; name) +``` + +Source that outputs a constant signal of `1` + +# Connectors +- `d` + Output [`DigitalPin`](@ref) +""" function Reset(; name) @named d = DigitalPin() @@ -29,6 +64,17 @@ function Reset(; name) ODESystem(eqs, t, [], [], systems=[d], name=name) end +""" +```julia +function Pulse(; name, duty_cycle=0.5, T=1.0) +``` + +Pulse output with specified `duty_cycle` and time period (`T`) + +# Connectors +- `d` + Output [`DigitalPin`](@ref) +""" function Pulse(; name, duty_cycle=0.5, T=1.0) @named d = DigitalPin() diff --git a/src/Electrical/Digital/tables.jl b/src/Electrical/Digital/tables.jl index 5f1999001..fe443b719 100644 --- a/src/Electrical/Digital/tables.jl +++ b/src/Electrical/Digital/tables.jl @@ -28,7 +28,7 @@ function _not(x) typeof(i) != Int && return X NotTable[i] end -@register _not(x) +@register_symbolic _not(x) # MISO AND gate AndTable = OffsetArray([ @@ -56,11 +56,11 @@ function _and(x...) end return y[end] end -@register _and(x...) -@register _and(a, b) -@register _and(a, b, c) -@register _and(a, b, c, d) -@register _and(a, b, c, d, e) +@register_symbolic _and(x...) +@register_symbolic _and(a, b) +@register_symbolic _and(a, b, c) +@register_symbolic _and(a, b, c, d) +@register_symbolic _and(a, b, c, d, e) # MISO OR gate OrTable = OffsetArray([ @@ -89,12 +89,12 @@ function _or(x...) end return y[end] end -@register _or(x...) -@register _or(a, b) -@register _or(a, b, c) -@register _or(a, b, c, d) -@register _or(a, b, c, d, e) -@register _or(a, b, c, d, e, f, g, h) +@register_symbolic _or(x...) +@register_symbolic _or(a, b) +@register_symbolic _or(a, b, c) +@register_symbolic _or(a, b, c, d) +@register_symbolic _or(a, b, c, d, e) +@register_symbolic _or(a, b, c, d, e, f, g, h) # MISO :XOR gate @@ -124,10 +124,10 @@ function _xor(x...) end return y[end] end -@register _xor(x...) -@register _xor(a, b) -@register _xor(a, b, c) -@register _xor(a, b, c, d) -@register _xor(a, b, c, d, e) +@register_symbolic _xor(x...) +@register_symbolic _xor(a, b) +@register_symbolic _xor(a, b, c) +@register_symbolic _xor(a, b, c, d) +@register_symbolic _xor(a, b, c, d, e) # TODO: revisit y[1] for all miso gates for 9-level logic \ No newline at end of file diff --git a/src/Electrical/Electrical.jl b/src/Electrical/Electrical.jl index 64a5d22e5..997205563 100644 --- a/src/Electrical/Electrical.jl +++ b/src/Electrical/Electrical.jl @@ -1,3 +1,7 @@ +""" +Library of electrical models. +This library contains electrical components to build up analog circuits. +""" module Electrical using ModelingToolkit, Symbolics, IfElse, OrdinaryDiffEq @@ -10,12 +14,19 @@ include("utils.jl") include("Analog/ideal_components.jl") include("Analog/sensors.jl") include("Analog/sources.jl") -include("Digital/components.jl") -include("Digital/gates.jl") -include("Digital/tables.jl") -include("Digital/sources.jl") +# include("Digital/components.jl") +# include("Digital/gates.jl") +# include("Digital/tables.jl") +# include("Digital/sources.jl") -export # Analog Components +# TODO: +# - digital +# - machines +# - multi-phase + +export #Interface + Pin, + # Analog Components Capacitor, Ground, Inductor, Resistor, Short, IdealOpAmp, # Analog Sensors @@ -24,17 +35,17 @@ export # Analog Components #Analog Sources ConstantVoltage, SineVoltage, StepVoltage, RampVoltage, SquareVoltage, TriangularVoltage, - CosineVoltage, DampedSineVoltage, + CosineVoltage, ExpSineVoltage, ConstantCurrent, SineCurrent, StepCurrent, RampCurrent, SquareCurrent, TriangularCurrent, - CosineCurrent, DampedSineCurrent, - connect, Pin + CosineCurrent, ExpSineCurrent + - # Digital Gates - And, Or, Not, Xor, Nand, Nor, Xnor, - # Digital components - HalfAdder, FullAdder, MUX, DEMUX, Encoder, Decoder, - # Digital Sources - DigitalPin, Pulse, PulseDiff + # # Digital Gates + # And, Or, Not, Xor, Nand, Nor, Xnor, + # # Digital components + # HalfAdder, FullAdder, MUX, DEMUX, Encoder, Decoder, + # # Digital Sources + # DigitalPin, Pulse, PulseDiff end diff --git a/src/Electrical/utils.jl b/src/Electrical/utils.jl index 9fa354e55..226be3cfc 100644 --- a/src/Electrical/utils.jl +++ b/src/Electrical/utils.jl @@ -1,64 +1,67 @@ @connector function Pin(;name) - @variables v(t) i(t) - ODESystem(Equation[], t, [v, i], [], name=name, defaults=Dict(v=>1.0, i=>1.0)) + sts = @variables begin + v(t) # Potential at the pin [V] + i(t), [connect=Flow] # Current flowing into the pin [A] + end + ODESystem(Equation[], t, sts, [], name=name, defaults=Dict(v=>1.0, i=>1.0)) end +@doc """ +```julia +@connector function Pin(; name) +``` -@connector function DigitalPin(; name) - @variables val(t) v(t) i(t) +A pin in an analog circuit. + +# States +- `v(t)`: [`V`] The voltage at this pin +- `i(t)`: [`A`] The current passing through this pin +""" Pin + +""" + OnePort(;name, v_start=0.0, i_start=0.0) + +Component with two electrical pins `p` and `n` and current `i` from `p` to `n`. + +# Parameters: +- `v_start`: [`V`] Initial voltage across the component +- `i_start`: [`A`] Initial current through the component +""" +function OnePort(;name, v_start=0.0, i_start=0.0) + @named p = Pin() + @named n = Pin() + sts = @variables begin + v(t)=v_start + i(t)=i_start + end eqs = [ - val ~ IfElse.ifelse((0.0 <= v) & (v <= 0.8) | (2.0 <= v) & (v <= 5.0), - IfElse.ifelse(v > 2.0, 1, 0), X) + v ~ p.v - n.v + 0 ~ p.i + n.i + i ~ p.i ] - ODESystem(Equation[], t, [val, v, i], [], defaults=Dict(val=>0, i=>0), name=name) + return compose(ODESystem(eqs, t, sts, []; name=name), p, n) end -abstract type ElectricalPin end -ModelingToolkit.promote_connect_rule(::Type{DigitalPin}, ::Type{Pin}) = ElectricalPin -ModelingToolkit.promote_connect_rule(::Type{Pin}, ::Type{DigitalPin}) = ElectricalPin -ModelingToolkit.promote_connect_rule(::Type{ElectricalPin}, ::Type{DigitalPin}) = ElectricalPin -ModelingToolkit.promote_connect_rule(::Type{ElectricalPin}, ::Type{Pin}) = ElectricalPin -function ModelingToolkit.connect(::Type{<:Pin}, ps...) - eqs = [ - 0 ~ sum(p->p.i, ps) # KCL - ] - # KVL - for i in 1:length(ps)-1 - push!(eqs, ps[i].v ~ ps[i+1].v) - end - return eqs -end - -function ModelingToolkit.connect(::Type{DigitalPin}, ps...) +@connector function DigitalPin(; name) + @variables val(t) v(t) i(t) eqs = [ - 0 ~ sum(p->p.i, ps) # KCL - ] - # KVL - for i in 1:length(ps)-1 - push!(eqs, ps[i].val ~ ps[i+1].val) - end - for i in 1:length(ps)-1 - push!(eqs, ps[i].v ~ ps[i+1].v) - end - return eqs + val ~ IfElse.ifelse((0.0 <= v) & (v <= 0.8) | (2.0 <= v) & (v <= 5.0), + IfElse.ifelse(v > 2.0, 1, 0), X) + ] + ODESystem(Equation[], t, [val, v, i], [], defaults = Dict(val => 0, i => 0), name = name) end -function ModelingToolkit.connect(::Type{ElectricalPin}, ps...) - eqs = [ - 0 ~ sum(p->p.i, ps) # KCL - ] +@doc """ +```julia +@connector function DigitalPin(; name) +``` - # KVL - digpins = ModelingToolkit.ODESystem[] - for p in ps - ModelingToolkit.get_connection_type(p) == DigitalPin && push!(digpins, p) - end - for i in 1:length(digpins)-1 - push!(eqs, digpins[i].val ~ digpins[i+1].val) - end - for i in 1:length(ps)-1 - push!(eqs, ps[i].v ~ ps[i+1].v) - end - return eqs -end +A pin in a digital circuit. + +# States +- `v(t)`: [`V`] The voltage at this pin +- `i(t)`: [`A`] The current passing through this pin +- `val(t)`: The binary value of the pin at this point. A voltage from `0V` to `0.8V` is a binary value + of `0`. A voltage in the range `2.0V` to `5.0V` is `1`. Any other value is `X`. +""" DigitalPin diff --git a/src/Magnetic/FluxTubes/FluxTubes.jl b/src/Magnetic/FluxTubes/FluxTubes.jl new file mode 100644 index 000000000..1657e499f --- /dev/null +++ b/src/Magnetic/FluxTubes/FluxTubes.jl @@ -0,0 +1,17 @@ +module FluxTubes +using ModelingToolkit +using ...Electrical: Pin + +@parameters t +D = Differential(t) + +export PositiveMagneticPort, NegativeMagneticPort, TwoPort +include("utils.jl") + +export Ground, Idle, Short, Crossing, ConstantPermeance, ConstantReluctance, EddyCurrent, ElectroMagneticConverter +include("basic.jl") + +export ConstantMagneticPotentialDifference, ConstantMagneticFlux +include("sources.jl") + +end #module \ No newline at end of file diff --git a/src/Magnetic/FluxTubes/basic.jl b/src/Magnetic/FluxTubes/basic.jl new file mode 100644 index 000000000..54671dd56 --- /dev/null +++ b/src/Magnetic/FluxTubes/basic.jl @@ -0,0 +1,149 @@ +""" + Ground(;name) + +Zero magnetic potential. +""" +function Ground(;name) + @named port = PositiveMagneticPort() + eqs = [port.V_m ~ 0] + ODESystem(eqs, t, [], [], systems=[port], name=name) +end + +""" + Idle(;name) + +Idle running branch. +""" +function Idle(;name) + @named two_port = TwoPort() + @unpack Phi = two_port + eqs = [ + Phi ~ 0, + ] + extend(ODESystem(eqs, t, [], [], systems=[], name=name), two_port) +end + +""" + Short(;name) + +Short cut branch. +""" +function Short(;name) + @named two_port = TwoPort() + @unpack V_m = two_port + eqs = [ + V_m ~ 0, + ] + extend(ODESystem(eqs, t, [], [], systems=[], name=name), two_port) +end + +""" + Crossing(;name) + +Crossing of two branches. + +This is a simple crossing of two branches. The ports port_p1 and port_p2 are connected, as well as port_n1 and port_n2. +""" +function Crossing(;name) + @named port_p1 = PositiveMagneticPort() + @named port_p2 = PositiveMagneticPort() + @named port_n1 = NegativeMagneticPort() + @named port_n2 = NegativeMagneticPort() + eqs = [ + connect(port_p1, port_p2), + connect(port_n1, port_n2), + ] + ODESystem(eqs, t, [], [], systems=[port_p1, port_p2, port_n1, port_n2], name=name) +end + +""" + ConstantPermeance(;name, G_m=1.0) + +Constant permeance. + +# Parameters: +- `G_m`: [H] Magnetic permeance +""" +function ConstantPermeance(;name, G_m=1.0) + @named two_port = TwoPort() + @unpack V_m, Phi = two_port + @parameters G_m=G_m + eqs = [ + Phi ~ G_m * V_m, + ] + extend(ODESystem(eqs, t, [], [G_m], name=name), two_port) +end + +""" + ConstantReluctance(;name, R_m=1.0) + +Constant reluctance. + +# Parameters: +- `R_m`: [H^-1] Magnetic reluctance +""" +function ConstantReluctance(;name, R_m=1.0) + @named two_port = TwoPort() + @unpack V_m, Phi = two_port + @parameters R_m=R_m + eqs = [ + V_m ~ Phi * R_m, + ] + extend(ODESystem(eqs, t, [], [R_m], name=name), two_port) +end + +""" + ElectroMagneticConverter(;name, N, Phi_start=0.0) + +Ideal electromagnetic energy conversion. + +The electromagnetic energy conversion is given by Ampere's law and Faraday's law respectively +V_m = N * i +N * dΦ/dt = -v + +# Parameters: +- `N`: Number of turns +- `Phi_start`: [Wb] Initial magnetic flux flowing into the port_p +""" +function ElectroMagneticConverter(;name, N, Phi_start=0.0) + @named port_p = PositiveMagneticPort() + @named port_n = NegativeMagneticPort() + @named p = Pin() + @named n = Pin() + + sts = @variables v(t) i(t) V_m(t) Phi(t)=Phi_start + pars = @parameters N=N + eqs = [ + v ~ p.v - n.v + 0 ~ p.i + n.i + i ~ p.i + V_m ~ port_p.V_m - port_n.V_m + 0 ~ port_p.Phi + port_n.Phi + Phi ~ port_p.Phi + #converter equations: + V_m ~ i * N # Ampere's law + D(Phi) ~ -v / N # Faraday's law + ] + ODESystem(eqs, t, sts, pars, systems=[port_p, port_n, p, n], name=name) +end + +""" + EddyCurrent(;name, rho=0.098e-6, l=1, A=1, Phi_start=0.0) + +For modelling of eddy current in a conductive magnetic flux tube. + +# Parameters: +- `rho`: [ohm * m] Resistivity of flux tube material (default: Iron at 20degC) +- `l`: [m] Average length of eddy current path +- `A`: [m^2] Cross sectional area of eddy current path +- `Phi_start`: [Wb] Initial magnetic flux flowing into the port_p +""" +function EddyCurrent(;name, rho=0.098e-6, l=1, A=1, Phi_start=0.0) + @named two_port = TwoPort(Phi_start=Phi_start) + @unpack V_m, Phi = two_port + @parameters R = rho * l / A # Electrical resistance of eddy current path + eqs = [ + D(Phi) ~ V_m * R, + ] + extend(ODESystem(eqs, t, [], [R], name=name), two_port) +end \ No newline at end of file diff --git a/src/Magnetic/FluxTubes/sources.jl b/src/Magnetic/FluxTubes/sources.jl new file mode 100644 index 000000000..1af781872 --- /dev/null +++ b/src/Magnetic/FluxTubes/sources.jl @@ -0,0 +1,37 @@ +""" +Constant magnetomotive force. + +Parameters: +- `V_m`: [A] Magnetic potential difference +""" +function ConstantMagneticPotentialDifference(;name, V_m=1.0) + port_p = PositiveMagneticPort() + port_n = NegativeMagneticPort() + @parameters V_m=V_m + @variables Phi(t) + eqs = [ + V_m ~ port_p.V_m - port_n.V_m + Phi ~ port_p.Phi + 0 ~ port_p.Phi + port_n.Phi + ] + ODESystem(eqs, t, [Phi], [V_m], systems=[port_p, port_n], name=name) +end + +""" +Source of constant magnetic flux. + +Parameters: +- `Phi`: [Wb] Magnetic flux +""" +function ConstantMagneticFlux(;name, Phi=1.0) + port_p = PositiveMagneticPort() + port_n = NegativeMagneticPort() + @parameters Phi=Phi + @variables V_m(t) + eqs = [ + V_m ~ port_p.V_m - port_n.V_m + Phi ~ port_p.Phi + 0 ~ port_p.Phi + port_n.Phi + ] + ODESystem(eqs, t, [Phi], [V_m], systems=[port_p, port_n], name=name) +end diff --git a/src/Magnetic/FluxTubes/utils.jl b/src/Magnetic/FluxTubes/utils.jl new file mode 100644 index 000000000..7b3d0b82c --- /dev/null +++ b/src/Magnetic/FluxTubes/utils.jl @@ -0,0 +1,37 @@ +@connector function MagneticPort(;name, V_m_start=0.0, Phi_start=0.0) + @variables V_m(t)=V_m_start # [Wb] Magnetic potential at the port + @variables Phi(t)=Phi_start [connect=Flow] # [A] Magnetic flux flowing into the port" + ODESystem(Equation[], t, [V_m, Phi], []; name=name) +end +Base.@doc "Port for a Magnetic system." MagneticPort + +""" +Positive magnetic port +""" +const PositiveMagneticPort = MagneticPort + +""" +Negative magnetic port +""" +const NegativeMagneticPort = MagneticPort + +""" + TwoPort(;name, V_m_start=0.0, Phi_start=0.0) + +Partial component with magnetic potential difference between two magnetic ports p and n and magnetic flux Phi from p to n. + +# Parameters: +- `V_m_start`: Initial magnetic potential difference between both ports +- `Phi_start`: Initial magnetic flux from port_p to port_n +""" +function TwoPort(;name, V_m_start=0.0, Phi_start=0.0) + @named port_p = PositiveMagneticPort() + @named port_n = NegativeMagneticPort() + @variables V_m(t)=V_m_start Phi(t)=Phi_start + eqs = [ + V_m ~ port_p.V_m - port_n.V_m + Phi ~ port_p.Phi + 0 ~ port_p.Phi + port_n.Phi + ] + ODESystem(eqs, t, [V_m, Phi], [], systems=[port_p, port_n], name=name) +end \ No newline at end of file diff --git a/src/Magnetic/Magnetic.jl b/src/Magnetic/Magnetic.jl index 87ef10efe..cd0295799 100644 --- a/src/Magnetic/Magnetic.jl +++ b/src/Magnetic/Magnetic.jl @@ -1,34 +1,12 @@ module Magnetic -using ModelingToolkit, Symbolics, IfElse, OrdinaryDiffEq +using ModelingToolkit -@parameters t -D = Differential(t) +# FluxTubes +include("FluxTubes/FluxTubes.jl") -include("utils.jl") +# QuasiStatic -include("FluxTubes/sensors.jl") - -include("QuasiStatic/FluxTubes/basic.jl") -include("QuasiStatic/FluxTubes/sources.jl") -include("QuasiStatic/FluxTubes/sensors.jl") - -export MagneticPort, - PositiveMagneticPort, - NegativeMagneticPort, - Ground, - TwoPortElementary, - TwoPortExtended, - TwoPort, - Idle, - Short, - Crossing, - ConstantPermeance, - ConstantReluctance, - ConstantMagneticPotentialDifference, - ConstantMagneticFlux, - AbsoluteSensor, - # FluxTubes sensors - MagneticFluxSensor, MagneticPotentialDifferenceSensor +# FundamentalWave end #module \ No newline at end of file diff --git a/src/Magnetic/QuasiStatic/FluxTubes/basic.jl b/src/Magnetic/QuasiStatic/FluxTubes/basic.jl deleted file mode 100644 index 79bb04380..000000000 --- a/src/Magnetic/QuasiStatic/FluxTubes/basic.jl +++ /dev/null @@ -1,57 +0,0 @@ -function Ground(;name) - @named port = PositiveMagneticPort() - eqs = [port.V_m ~ 0] - ODESystem(eqs, t, [], [], systems=[port], name=name) -end - -function Idle(;name) - @named two_port = TwoPort() - @unpack Phi = two_port - eqs = [ - Phi ~ 0, - ] - extend(ODESystem(eqs, t, [Phi], [], systems=[], name=name), two_port) -end - -function Short(;name) - port_p = PositiveMagneticPort() - port_n = NegativeMagneticPort() - eqs = [ - connect(port_p, port_n), - ] - ODESystem(eqs, t, [], [], systems=[port_p, port_n], name=name) -end - -function Crossing(;name) - port_p1 = PositiveMagneticPort() - port_p2 = PositiveMagneticPort() - port_n1 = NegativeMagneticPort() - port_n2 = NegativeMagneticPort() - eqs = [ - connect(port_p1, port_p2), - connect(port_n1, port_n2), - ] - ODESystem(eqs, t, [], [], systems=[port_p1, port_p2, port_n1, port_n2], name=name) -end - -function ConstantPermeance(;name, G_m=1.0) - val = G_m - @named two_port = TwoPort() - @unpack V_m, Phi = two_port - @variables G_m(t) - eqs = [ - Phi ~ G_m * V_m, - ] - extend(ODESystem(eqs, t, [V_m, Phi, G_m], [], systems=[], defaults=Dict(G_m => val), name=name), two_port) -end - -function ConstantReluctance(;name, R_m=1.0) - val = R_m - @named two_port = TwoPort() - @unpack V_m, Phi = two_port - @variables R_m(t) - eqs = [ - V_m ~ Phi * R_m, - ] - extend(ODESystem(eqs, t, [V_m, Phi, R_m], [], systems=[], defaults=Dict(R_m => val), name=name), two_port) -end diff --git a/src/Magnetic/QuasiStatic/FluxTubes/sources.jl b/src/Magnetic/QuasiStatic/FluxTubes/sources.jl deleted file mode 100644 index 33a19865f..000000000 --- a/src/Magnetic/QuasiStatic/FluxTubes/sources.jl +++ /dev/null @@ -1,27 +0,0 @@ -function ConstantMagneticPotentialDifference(;name, V_m=1.0) - val = V_m - @named two_port_elementary = TwoPortElementary() - @unpack port_p, port_n = two_port_elementary - @parameters V_m - @variables Phi(t) - eqs = [ - V_m ~ port_p.V_m - port_n.V_m, - Phi ~ port_p.Phi, - 0 ~ port_p.Phi + port_n.Phi, - ] - extend(ODESystem(eqs, t, [Phi], [V_m], systems=[port_p, port_n], defaults=Dict(V_m => val), name=name), two_port_elementary) -end - -function ConstantMagneticFlux(;name, Phi=1.0) - val = Phi - @named two_port_elementary = TwoPortElementary() - @unpack port_p, port_n = two_port_elementary - @parameters Phi - @variables V_m(t) - eqs = [ - V_m ~ port_p.V_m - port_n.V_m, - Phi ~ port_p.Phi, - 0 ~ port_p.Phi + port_n.Phi, - ] - extend(ODESystem(eqs, t, [V_m], [Phi], systems=[port_p, port_n], defaults=Dict(Phi => val), name=name), two_port_elementary) -end diff --git a/src/Magnetic/utils.jl b/src/Magnetic/utils.jl deleted file mode 100644 index 0bf2ba7ae..000000000 --- a/src/Magnetic/utils.jl +++ /dev/null @@ -1,59 +0,0 @@ -@connector function MagneticPort(;name, complex=false) - if complex - V_m, Phi = @variables V_m(t)::Complex Phi(t)::Complex - else - V_m, Phi = @variables V_m(t) Phi(t) - end - ODESystem(Equation[], t, [V_m, Phi], [], name=name, defaults=Dict(V_m=>0.0, Phi=>0.0)) -end - -function ModelingToolkit.connect(::Type{<:MagneticPort}, ps...) - eqs = [ - 0 ~ sum(p->p.Phi, ps) # Gauss's law for magnetism - ] - - for i in 1:length(ps)-1 - push!(eqs, ps[i].V_m ~ ps[i+1].V_m) - end - - return eqs -end - -const PositiveMagneticPort = MagneticPort -const NegativeMagneticPort = MagneticPort - -function TwoPortElementary(;name, complex=false) - @named port_p = PositiveMagneticPort(;complex=complex) - @named port_n = NegativeMagneticPort(;complex=complex) - ODESystem(Equation[], t, [], [], systems=[port_p, port_n], name=name) -end - -function TwoPortExtended(;name, complex=false) - @named two_port_elementary = TwoPortElementary(complex=complex) - @unpack port_p, port_n = two_port_elementary - @variables V_m(t) Phi(t) - eqs = [ - V_m ~ port_p.V_m - port_n.V_m, - Phi ~ port_p.Phi, - ] - extend(ODESystem(eqs, t, [V_m, Phi], [], systems=[port_p, port_n], name=name), two_port_elementary) -end - -function TwoPort(;name, complex=false) - @named two_port_extended = TwoPortExtended(;complex=complex) - @unpack port_p, port_n = two_port_extended - eqs = [ - 0 ~ port_p.Phi + port_n.Phi, - ] - extend(ODESystem(eqs, t, [], [], systems=[port_p, port_n], name=name), two_port_extended) -end - -function AbsoluteSensor(;name) - @variables omega - @named port = PositiveMagneticPort(;complex=true) - eqs = [ - D(port.reference.gamma) ~ omega, - port.Phi ~ Complex(0); - ] - ODESystem(eqs, t, [omega, port.Phi, port.reference.gamma], [], systems=[], name=name) -end \ No newline at end of file diff --git a/src/Mechanical/Mechanical.jl b/src/Mechanical/Mechanical.jl new file mode 100644 index 000000000..3caaf5ae7 --- /dev/null +++ b/src/Mechanical/Mechanical.jl @@ -0,0 +1,11 @@ +""" +Library of mechanical models. +""" +module Mechanical + +using ModelingToolkit + +include("Rotational/Rotational.jl") + +end + diff --git a/src/Mechanical/Rotational/Rotational.jl b/src/Mechanical/Rotational/Rotational.jl new file mode 100644 index 000000000..0a4dd4798 --- /dev/null +++ b/src/Mechanical/Rotational/Rotational.jl @@ -0,0 +1,22 @@ +""" +Library to model 1-dimensional, rotational mechanical systems +""" +module Rotational + +using ModelingToolkit, Symbolics, IfElse, OrdinaryDiffEq +using OffsetArrays +using ...Blocks: RealInput, RealOutput + +@parameters t +D = Differential(t) + +export Flange +include("utils.jl") + +export Fixed, Inertia, Spring, Damper, IdealGear +include("components.jl") + +export Torque +include("sources.jl") + +end \ No newline at end of file diff --git a/src/Mechanical/Rotational/components.jl b/src/Mechanical/Rotational/components.jl new file mode 100644 index 000000000..d865eca36 --- /dev/null +++ b/src/Mechanical/Rotational/components.jl @@ -0,0 +1,110 @@ +""" + Fixed(;name, phi0=0.0) + +Flange fixed in housing at a given angle. + +# Parameters: +- `phi0`: [rad] Fixed offset angle of housing +""" +function Fixed(;name, phi0=0.0) + @named flange = Flange() + @parameters phi0=phi0 + eqs = [flange.phi ~ phi0] + return compose(ODESystem(eqs, t, [], [phi0]; name=name), flange) +end + +""" + Inertia(;name, J=1.0, phi_start=0.0, w_start=0.0, a_start=0.0) + +1D-rotational component with inertia. + +# Parameters: +- `J`: [kg·m²] Moment of inertia +- `phi_start`: [rad] Initial value of absolute rotation angle of component +- `w_start`: [rad/s] Initial value of absolute angular velocity of component +- `a_start`: [rad/s²] Initial value of absolute angular acceleration of component + +# States: +- `phi`: [rad] Absolute rotation angle of component +- `w`: [rad/s] Absolute angular velocity of component (= der(phi)) +- `a`: [rad/s²] Absolute angular acceleration of component (= der(w)) +""" +function Inertia(;name, J=1.0, phi_start=0.0, w_start=0.0, a_start=0.0) + @named flange_a = Flange() + @named flange_b = Flange() + @parameters J=J + sts = @variables begin + phi(t)=phi_start + w(t)=w_start + a(t)=a_start + end + eqs = [ + phi ~ flange_a.phi + phi ~ flange_b.phi + D(phi) ~ w + D(w) ~ a + J*a ~ flange_a.tau + flange_b.tau + ] + return compose(ODESystem(eqs, t, sts, [J]; name=name), flange_a, flange_b) +end + +""" + Spring(;name, c, phi_rel0=0.0) + +Linear 1D rotational spring + +# Parameters: +- `c`: [N.m/rad] Spring constant +- `phi_rel0`: Unstretched spring angle +""" +function Spring(;name, c=1.0e5, phi_rel0=0.0) + @named partial_comp = PartialCompliant() + @unpack phi_rel, tau = partial_comp + pars = @parameters begin + c=c + phi_rel0=phi_rel0 + end + eqs = [tau ~ c*(phi_rel - phi_rel0)] + extend(ODESystem(eqs, t, [], pars; name=name), partial_comp) +end + +""" + Damper(;name, d=0.0) + +Linear 1D rotational damper + +# Parameters: +- `d`: [N.m.s/rad] Damping constant +""" +function Damper(;name, d=0.0) + @named partial_comp = PartialCompliantWithRelativeStates() + @unpack w_rel, tau = partial_comp + pars = @parameters d=d + eqs = [tau ~ d*w_rel] + extend(ODESystem(eqs, t, [], pars; name=name), partial_comp) +end + +""" + IdealGear(;name, ratio, use_support=false) + +Ideal gear without inertia. + +This element characterizes any type of gear box which is fixed in the ground and which has one driving shaft and one driven shaft. + +# Parameters: +- `ratio`: Transmission ratio (flange_a.phi/flange_b.phi) +- `use_support`: If support flange enabled, otherwise implicitly grounded +""" +function IdealGear(;name, ratio, use_support=false) + @named partial_element = PartialElementaryTwoFlangesAndSupport2(use_support=use_support) + @unpack phi_support, flange_a, flange_b = partial_element + @parameters ratio=ratio + sts = @variables phi_a(t)=0.0 phi_b(t)=0.0 + eqs = [ + phi_a ~ flange_a.phi - phi_support + phi_b ~ flange_b.phi - phi_support + phi_a ~ ratio*phi_b + 0 ~ ratio*flange_a.tau + flange_b.tau + ] + extend(ODESystem(eqs, t, sts, [ratio]; name=name), partial_element) +end \ No newline at end of file diff --git a/src/Mechanical/Rotational/sources.jl b/src/Mechanical/Rotational/sources.jl new file mode 100644 index 000000000..28e6f2503 --- /dev/null +++ b/src/Mechanical/Rotational/sources.jl @@ -0,0 +1,12 @@ +""" + Torque(;name) + +Input signal acting as external torque on a flange +""" +function Torque(;name, use_support=false) + @named partial_element = PartialElementaryOneFlangeAndSupport2(use_support=use_support) + @unpack flange = partial_element + @named tau = RealInput() # Accelerating torque acting at flange (= -flange.tau) + eqs = [flange.tau ~ -tau.u] + return extend(ODESystem(eqs, t, [], []; name=name, systems=[tau]), partial_element) +end \ No newline at end of file diff --git a/src/Mechanical/Rotational/utils.jl b/src/Mechanical/Rotational/utils.jl new file mode 100644 index 000000000..730802edf --- /dev/null +++ b/src/Mechanical/Rotational/utils.jl @@ -0,0 +1,154 @@ +@connector function Flange(;name) + sts = @variables begin + phi(t) + tau(t), [connect=Flow] + end + ODESystem(Equation[], t, sts, [], name=name, defaults=Dict(phi=>0.0, tau=>0.0)) +end +Base.@doc """ + Flange(;name) + +1-dim. rotational flange of a shaft. + +# States: +- `phi`: [rad] Absolute rotation angle of flange +- `tau`: [N.m] Cut torque in the flange +""" Flange + +@connector function Support(;name) + sts = @variables begin + phi(t) + tau(t), [connect=Flow] + end + ODESystem(Equation[], t, sts, [], name=name, defaults=Dict(phi=>0.0, tau=>0.0)) +end +Base.@doc """ + Support(;name) + +Support/housing of a 1-dim. rotational shaft + +# States: +- `phi`: [rad] Absolute rotation angle of the support/housing +- `tau`: [N.m] Cut torque in the support/housing +""" Support + +""" + PartialCompliant(;name, phi_rel_start=0.0, tau_start=0.0) + +Partial model for the compliant connection of two rotational 1-dim. shaft flanges. + +# Parameters: +- `phi_rel_start`: [rad] Initial relative rotation angle +- `tau_start`: [N.m] Initial torque between flanges + +# States: +- `phi_rel`: [rad] Relative rotation angle (= flange_b.phi - flange_a.phi) +- `tau`: [N.m] Torque between flanges (= flange_b.tau) +""" +function PartialCompliant(;name, phi_rel_start=0.0, tau_start=0.0) + @named flange_a = Flange() + @named flange_b = Flange() + sts = @variables begin + phi_rel(t)=phi_rel_start + tau(t)=tau_start + end + eqs = [ + phi_rel ~ flange_b.phi - flange_a.phi + flange_b.tau ~ tau + flange_a.tau ~ -tau + ] + return compose(ODESystem(eqs, t, sts, []; name=name), flange_a, flange_b) +end + +""" + PartialCompliantWithRelativeStates(;name, phi_rel_start=0.0, tau_start=0.0) + +Partial model for the compliant connection of two rotational 1-dim. shaft flanges where the relative angle and speed are used as preferred states + +# Parameters: +- `phi_rel_start`: [rad] Initial relative rotation angle +- `w_rel_start`: [rad/s] Initial relative angular velocity (= der(phi_rel)) +- `a_rel_start`: [rad/s²] Initial relative angular acceleration (= der(w_rel)) +- `tau_start`: [N.m] Initial torque between flanges + +# States: +- `phi_rel`: [rad] Relative rotation angle (= flange_b.phi - flange_a.phi) +- `w_rel`: [rad/s] Relative angular velocity (= der(phi_rel)) +- `a_rel`: [rad/s²] Relative angular acceleration (= der(w_rel)) +- `tau`: [N.m] Torque between flanges (= flange_b.tau) +""" +function PartialCompliantWithRelativeStates(;name, phi_rel_start=0.0, w_start=0.0, a_start=0.0, tau_start=0.0) + @named flange_a = Flange() + @named flange_b = Flange() + sts = @variables begin + phi_rel(t)=phi_rel_start + w_rel(t)=w_start + a_rel(t)=a_start + tau(t)=tau_start + end + eqs = [ + phi_rel ~ flange_b.phi - flange_a.phi + D(phi_rel) ~ w_rel + D(w_rel) ~ a_rel + flange_b.tau ~ tau + flange_a.tau ~ -tau + ] + return compose(ODESystem(eqs, t, sts, []; name=name), flange_a, flange_b) +end + +""" + PartialElementaryOneFlangeAndSupport2(;name, use_support=false) + +Partial model for a component with one rotational 1-dim. shaft flange and a support used for textual modeling, i.e., for elementary models + +# Parameters: +- `use_support`: If support flange enabled, otherwise implicitly grounded + +# States: +- `phi_support`: [rad] Absolute angle of support flange" +""" +function PartialElementaryOneFlangeAndSupport2(;name, use_support=false) + @named flange = Flange() + sys = [flange] + @variables phi_support(t) + if use_support + @named support = Support() + eqs = [ + support.phi ~ phi_support + support.tau ~ -flange.tau + ] + push!(sys, support) + else + eqs = [phi_support ~ 0] + end + return compose(ODESystem(eqs, t, [phi_support], []; name=name), sys) +end + +""" + PartialElementaryTwoFlangesAndSupport2(;name, use_support=false) + +Partial model for a component with two rotational 1-dim. shaft flanges and a support used for textual modeling, i.e., for elementary models + +# Parameters: +- `use_support`: If support flange enabled, otherwise implicitly grounded + +# States: +- `phi_support`: [rad] Absolute angle of support flange" +""" +function PartialElementaryTwoFlangesAndSupport2(;name, use_support=false) + @named flange_a = Flange() + @named flange_b = Flange() + sys = [flange_a, flange_b] + @variables phi_support(t)=0.0 + if use_support + @named support = Support() + eqs = [ + support.phi ~ phi_support + support.tau ~ -flange_a.tau - flange_b.tau + ] + push!(sys, support) + else + eqs = [phi_support ~ 0] + end + return compose(ODESystem(eqs, t, [phi_support], []; name=name), sys) +end \ No newline at end of file diff --git a/src/ModelingToolkitStandardLibrary.jl b/src/ModelingToolkitStandardLibrary.jl index e5666f2e1..e8a24088d 100644 --- a/src/ModelingToolkitStandardLibrary.jl +++ b/src/ModelingToolkitStandardLibrary.jl @@ -1,8 +1,9 @@ module ModelingToolkitStandardLibrary include("Blocks/Blocks.jl") +include("Mechanical/Mechanical.jl") include("Electrical/Electrical.jl") -#include("Magnetic/Magnetic.jl") +include("Magnetic/Magnetic.jl") include("Thermal/Thermal.jl") end diff --git a/src/Thermal/HeatTransfer/ideal_components.jl b/src/Thermal/HeatTransfer/ideal_components.jl index 234485e8f..36f373671 100644 --- a/src/Thermal/HeatTransfer/ideal_components.jl +++ b/src/Thermal/HeatTransfer/ideal_components.jl @@ -1,121 +1,162 @@ -function ThermalGround(; name) - @named hp = HeatPort() - eqs = [hp.T ~ 0] - ODESystem(eqs, t, systems=[hp], name=name) -end - -function HeatCapacitor(; name, C=1.0) - c_th = C - - @named hp = HeatPort() - @parameters C - @variables T(t) dt(t) +""" + HeatCapacitor(; name, C=1.0, T_start=293.15 + 20) + +Lumped thermal element storing heat + +# Parameters: +- `C`: [J/K] Heat capacity of element (= cp*m) +- `T_start`: Initial temperature of element + +# States: +- `T`: [K] Temperature of element +- `der_T`: [K/s] Time derivative of temperature +""" +function HeatCapacitor(; name, C=1.0, T_start=293.15 + 20) + @named port = HeatPort() + @parameters C=C + sts = @variables begin + T(t)=T_start + der_T(t) + end D = Differential(t) eqs = [ - T ~ hp.T - dt ~ D(T) - D(T) ~ hp.Q_flow / C - ] - ODESystem(eqs, t, [T, dt], [C], systems=[hp], defaults=Dict(C => c_th), name=name) + T ~ port.T + der_T ~ D(T) + D(T) ~ port.Q_flow / C + ] + ODESystem(eqs, t, sts, [C]; systems=[port], name=name) end -function ThermalConductor(; name, G=1.0) - g_th = G +""" + ThermalConductor(;name, G=1.0) - @named hp1 = HeatPort() - @named hp2 = HeatPort() - @parameters G - @variables Q_flow(t) T(t) +Lumped thermal element transporting heat without storing it. +# Parameters: +- `G`: [W/K] Constant thermal conductance of material +""" +function ThermalConductor(;name, G=1.0) + @named element1d = Element1D() + @unpack Q_flow, dT = element1d + pars = @parameters G=G eqs = [ - T ~ hp1.T - hp2.T - Q_flow ~ G*T - Q_flow ~ hp1.Q_flow - -Q_flow ~ hp2.Q_flow + Q_flow ~ G * dT ] - ODESystem(eqs, t, [Q_flow, T], [G], systems=[hp1, hp2], defaults=Dict(G => g_th), name=name) + extend(ODESystem(eqs, t, [], pars; name=name), element1d) end -function ThermalResistor(; name, R=1.0) - r_th = R - @named hp1 = HeatPort() - @named hp2 = HeatPort() - @parameters R - @variables Q_flow(t) T(t) +""" + ThermalResistor(; name, R=1.0) + +Lumped thermal element transporting heat without storing it. +# Parameters: +- `R`: [K/W] Constant thermal resistance of material +""" +function ThermalResistor(; name, R=1.0) + @named element1d = Element1D() + @unpack Q_flow, dT = element1d + pars = @parameters R=R eqs = [ - T ~ R*Q_flow - T ~ hp1.T - hp2.T - hp1.Q_flow ~ Q_flow - hp2.Q_flow ~ -Q_flow + dT ~ R * Q_flow ] - ODESystem(eqs, t, [Q_flow, T], [R], systems=[hp1, hp2], defaults=Dict(R => r_th), name=name) + + extend(ODESystem(eqs, t, [], pars; name=name), element1d) end -function ConvectiveConductor(; name, G=1.0) - g_c = G +""" + ConvectiveConductor(; name, G=1.0) - @named solidport = HeatPort() - @named fluidport = HeatPort() - @parameters G # Convective thermal conductance - @variables Q_flow(t) dT(t) +Lumped thermal element for heat convection. +# Parameters: +- `G`: [W/K] Convective thermal conductance + +# States: +- `dT`: [K] Temperature difference across the component solid.T - fluid.T +- `Q_flow`: [W] Heat flow rate from solid -> fluid +""" +function ConvectiveConductor(; name, G=1.0) + @named solid = HeatPort() + @named fluid = HeatPort() + @parameters G=G + sts = @variables Q_flow(t) dT(t) eqs = [ - dT ~ solidport.T - fluidport.T - solidport.Q_flow ~ Q_flow - fluidport.Q_flow ~ -Q_flow + dT ~ solid.T - fluid.T + solid.Q_flow ~ Q_flow + fluid.Q_flow ~ -Q_flow dT ~ G*Q_flow ] - ODESystem(eqs, t, [Q_flow, dT], [G], systems=[solidport, fluidport], defaults=Dict(G => g_c), name=name) + ODESystem(eqs, t, sts, [G]; systems=[solid, fluid], name=name) end -function ConvectiveResistor(; name, R=1.0) - r_c = R +""" + ConvectiveResistor(; name, R=1.0) + +Lumped thermal element for heat convection. +# Parameters: +- `R`: [K/W] Constant thermal resistance of material + +# States: +- `dT`: [K] Temperature difference across the component solid.T - fluid.T +- `Q_flow`: [W] Heat flow rate from solid -> fluid +""" +function ConvectiveResistor(; name, R=1.0) @named solidport = HeatPort() @named fluidport = HeatPort() - @parameters R # Convective thermal resistance - @variables Q_flow(t) dT(t) - + @parameters R=R + sts = @variables Q_flow(t) dT(t) eqs = [ dT ~ solidport.T - fluidport.T solidport.Q_flow ~ Q_flow fluidport.Q_flow ~ -Q_flow dT ~ R*Q_flow ] - ODESystem(eqs, t, [Q_flow, dT], [R], systems=[solidport, fluidport], defaults=Dict(R => r_c), name=name) + ODESystem(eqs, t, sts, [R]; systems=[solidport, fluidport], name=name) end +""" + BodyRadiation(; name, G=1.0) + +Lumped thermal element for radiation heat transfer. + +# Parameters: +- `G`: [m^2] Net radiation conductance between two surfaces +""" function BodyRadiation(; name, G=1.0) - g_r = G - σ = 5.6703744191844294e-8 # Stefan-Boltzmann constant - - @named hp1 = HeatPort() - @named hp2 = HeatPort() - @parameters G # Net radiation conductance between two surfaces - @variables Q_flow(t) + sigma = 5.6703744191844294e-8 # Stefan-Boltzmann constant TODO: extract into physical constants module or use existing one + @named element1d = Element1D() + @unpack Q_flow, dT = element1d + @unpack port_a, port_b = element1d + pars = @parameters G=G eqs = [ - Q_flow ~ G*σ*(hp1.T^4 - hp2.T^4) + Q_flow ~ G * sigma * (port_a.T^4 - port_b.T^4) ] - ODESystem(eqs, t, [Q_flow], [G], systems=[hp1, hp2], defaults=Dict(G => g_r), name=name) + + extend(ODESystem(eqs, t, [], pars; name=name), element1d) end -function ThermalCollector(; name, N=1) - hp = [] - for i in 1:N - _hp = HeatPort(name=Symbol(:hp, i)) - push!(hp, _hp) - end - @named collector_port = HeatPort() +""" + ThermalCollector(; name, m=1) + +Collects m heat flows +This is a model to collect the heat flows from `m` heatports to one single heatport. +# Parameters: +- `m`: Number of heat ports (e.g. m=2: `port_a1`, `port_a2`) +""" +function ThermalCollector(; name, m=1) + port_a = [HeatPort(name=Symbol(:port_a, i)) for i in 1:m] + @named port_b = HeatPort() eqs = [ - collector_port.Q_flow + sum(k -> k.Q_flow, hp) ~ 0 - collector_port.T ~ hp[1].T + port_b.Q_flow + sum(k -> k.Q_flow, port_a) ~ 0 + port_b.T ~ port_a[1].T ] - for i in 1:N-1 - push!(eqs, hp[i].T ~ hp[i+1].T) + for i in 1:m-1 + push!(eqs, port_a[i].T ~ port_a[i+1].T) end - ODESystem(eqs, t, [], [], systems=[hp..., collector_port], name=name) + ODESystem(eqs, t, [], []; systems=[port_a..., port_b], name=name) end diff --git a/src/Thermal/HeatTransfer/sensors.jl b/src/Thermal/HeatTransfer/sensors.jl index 28948ddaa..d7a746e74 100644 --- a/src/Thermal/HeatTransfer/sensors.jl +++ b/src/Thermal/HeatTransfer/sensors.jl @@ -1,36 +1,60 @@ +""" + TemperatureSensor(; name) + +Absolute temperature sensor in kelvin. + +This is an ideal absolute temperature sensor which returns the temperature of the connected port in kelvin as an output +signal. The sensor itself has no thermal interaction with whatever it is connected to. Furthermore, no thermocouple-like +lags are associated with this sensor model. +""" function TemperatureSensor(; name) - @named hp = HeatPort() - @variables T(t) - + @named port = HeatPort() + @variables T(t) # [K] Absolute temperature eqs = [ - T ~ hp.T - hp.Q_flow ~ 0 + T ~ port.T + port.Q_flow ~ 0 ] - ODESystem(eqs, t, [T], [], systems=[hp], name=name) + ODESystem(eqs, t, [T], [], systems=[port], name=name) end +""" + RelativeTemperatureSensor(; name) + +Relative Temperature sensor. + +The relative temperature `port_a.T - port_b.T` is determined between the two ports of this component and is provided as +output signal in kelvin. +""" function RelativeTemperatureSensor(; name) - @named hp1 = HeatPort() - @named hp2 = HeatPort() - @variables T(t) - + @named port_a = HeatPort() + @named port_b = HeatPort() + @variables T(t) # [K] Relative temperature a.T - b.T eqs = [ - T ~ hp1.T - hp2.T - hp1.Q_flow ~ 0 - hp2.Q_flow ~ 0 + T ~ port_a.T - port_b.T + port_a.Q_flow ~ 0 + port_b.Q_flow ~ 0 ] - ODESystem(eqs, t, [T], [], systems=[hp1, hp2], name=name) + ODESystem(eqs, t, [T], [], systems=[port_a, port_b], name=name) end +""" + HeatFlowSensor(; name) + +Heat flow rate sensor. + +This model is capable of monitoring the heat flow rate flowing through this component. The sensed value of heat flow rate +is the amount that passes through this sensor while keeping the temperature drop across the sensor zero. This is an ideal +model so it does not absorb any energy and it has no direct effect on the thermal response of a system it is included in. +The output signal is positive, if the heat flows from `port_a` to `port_b`. +""" function HeatFlowSensor(; name) - @named hp1 = HeatPort() - @named hp2 = HeatPort() - @variables Q_flow(t) - + @named port_a = HeatPort() + @named port_b = HeatPort() + @variables Q_flow(t) # [W] Heat flow from port a to port b eqs = [ - hp1.T ~ hp2.T - hp1.Q_flow + hp2.Q_flow ~ 0 - Q_flow ~ hp1.Q_flow + port_a.T ~ port_b.T + port_a.Q_flow + port_b.Q_flow ~ 0 + Q_flow ~ port_a.Q_flow ] - ODESystem(eqs, t, [Q_flow], [], systems=[hp1, hp2], name=name) + ODESystem(eqs, t, [Q_flow], [], systems=[port_a, port_b], name=name) end diff --git a/src/Thermal/HeatTransfer/sources.jl b/src/Thermal/HeatTransfer/sources.jl index eb7fa92e5..be3297ecd 100644 --- a/src/Thermal/HeatTransfer/sources.jl +++ b/src/Thermal/HeatTransfer/sources.jl @@ -1,23 +1,46 @@ -function FixedHeatFlow(; name, Q_flow=1.0, T₀=293.15, α=0.0) - qflow, tem₀, alpha = Q_flow, T₀, α - - @parameters Q_flow T₀ α - @named hp = HeatPort() +""" + FixedHeatFlow(; name, Q_flow=1.0, T_ref=293.15, alpha=0.0) + +Fixed heat flow boundary condition. + +This model allows a specified amount of heat flow rate to be "injected" into a thermal system at a given port. +The constant amount of heat flow rate `Q_flow` is given as a parameter. The heat flows into the component to which +the component FixedHeatFlow is connected, if parameter `Q_flow` is positive. + +# Parameters: +- `Q_flow`: [W] Fixed heat flow rate at port +- `T_ref`: [K] Reference temperature +- `alpha`: [1/K] Temperature coefficient of heat flow rate +""" +function FixedHeatFlow(; name, Q_flow=1.0, T_ref=293.15, alpha=0.0) + pars = @parameters begin + Q_flow=Q_flow + T_ref=T_ref + alpha=alpha + end + @named port = HeatPort() eqs = [ - hp.Q_flow ~ -Q_flow * (1 + α*(hp.T - T₀)) + port.Q_flow ~ -Q_flow * (1 + alpha * (port.T - T_ref)) ] - ODESystem(eqs, t, [], [Q_flow, T₀, α], systems=[hp], defaults=Dict(zip((Q_flow, T₀, α), (qflow, tem₀, alpha))), name=name) + ODESystem(eqs, t, [], pars; systems=[port], name=name) end -function FixedTemperature(; name, T=0.0) - tem = T +""" + FixedTemperature(; name, T=0.0) + +Fixed temperature boundary condition in kelvin. - @named hp = HeatPort() - @parameters T +This model defines a fixed temperature T at its port in kelvin, i.e., it defines a fixed temperature as a boundary condition. +# Parameters: +- `T`: [K] Fixed temperature boundary condition +""" +function FixedTemperature(; name, T=0.0) + @named port = HeatPort() + pars = @parameters T=T eqs = [ - hp.T ~ T + port.T ~ T ] - ODESystem(eqs, t, [], [T], systems=[hp], defaults=Dict(T => tem), name=name) + ODESystem(eqs, t, [], pars; systems=[port], name=name) end diff --git a/src/Thermal/Thermal.jl b/src/Thermal/Thermal.jl index 0880f1234..5cde55cf0 100644 --- a/src/Thermal/Thermal.jl +++ b/src/Thermal/Thermal.jl @@ -1,22 +1,31 @@ +""" +Library of thermal system components to model heat transfer. +""" module Thermal using ModelingToolkit, Symbolics, IfElse, OrdinaryDiffEq @parameters t D = Differential(t) +include("utils.jl") + +# Library of 1-dimensional heat transfer with lumped elements include("HeatTransfer/ideal_components.jl") include("HeatTransfer/sensors.jl") include("HeatTransfer/sources.jl") -include("utils.jl") - -export # Thermal Components +# Simple components for 1-dimensional incompressible thermo-fluid flow models +# TODO: +# - FluidHeatFlow + +export # Interface + HeatPort, + # Thermal Components BodyRadiation, ConvectiveConductor, ConvectiveResistor, HeatCapacitor, ThermalConductor, ThermalResistor, ThermalCollector, # Thermal Sensors RelativeTemperatureSensor, HeatFlowSensor, TemperatureSensor, # Thermal Sources - FixedHeatFlow, FixedTemperature, ThermalGround, - connect, HeatPort + FixedHeatFlow, FixedTemperature, ThermalGround end \ No newline at end of file diff --git a/src/Thermal/utils.jl b/src/Thermal/utils.jl index dacdc256f..b355ab89c 100644 --- a/src/Thermal/utils.jl +++ b/src/Thermal/utils.jl @@ -1,15 +1,49 @@ -@connector function HeatPort(; name) - @variables T(t), Q_flow(t) # Temperature and Heat-flow-rate +@connector function HeatPort(; name, T_start=273.15 + 20.0, Q_flow_start=0.0) + @variables T(t)=T_start + @variables Q_flow(t)=Q_flow_start [connect = Flow] ODESystem(Equation[], t, [T, Q_flow], [], name=name) end +Base.@doc """ + HeatPort(; name, T_start=273.15 + 20.0, Q_flow_start=0.0) -function ModelingToolkit.connect(::Type{<:HeatPort}, ps...) - eqs = [ - 0 ~ sum(p->p.Q_flow, ps) - ] - for i in 1:length(ps)-1 - push!(eqs, ps[i].T ~ ps[i+1].T) - end +Port for a thermal system. - return eqs -end +# Parameters: +- `T_start`: [K] Temperature of the port +- `Q_flow_start`: [W] Heat flow rate at the port + +# States: +- `T`: [K] Temperature of the port +- `Q_flow`: [W] Heat flow rate at the port +""" HeatPort + +""" + Element1D(;name, dT0=0.0, Q_flow0=0.0) + +This partial model contains the basic connectors and variables to allow heat transfer models to be created that do not +store energy. This model defines and includes equations for the temperature drop across the element, `dT`, and the heat +flow rate through the element from `port_a` to `port_b`, `Q_flow`. + +# Parameters: +- `dT_start`: [K] Initial temperature difference across the component a.T - b.T +- `Q_flow_start`: [W] Initial heat flow rate from port a -> port b + +# States: +- `dT`: [K] Temperature difference across the component a.T - b.T +- `Q_flow`: [W] Heat flow rate from port a -> port b +""" +function Element1D(;name, dT_start=0.0, Q_flow_start=0.0) + @named port_a = HeatPort() + @named port_b = HeatPort() + sts = @variables begin + dT(t)=dT_start + Q_flow(t)=Q_flow_start + end + eqs = [ + dT ~ port_a.T - port_b.T + port_a.Q_flow ~ Q_flow + port_a.Q_flow + port_b.Q_flow ~ 0 + ] + + return compose(ODESystem(eqs, t, sts, []; name=name), port_a, port_b) +end \ No newline at end of file diff --git a/test/Blocks/continuous.jl b/test/Blocks/continuous.jl new file mode 100644 index 000000000..cadc5fb08 --- /dev/null +++ b/test/Blocks/continuous.jl @@ -0,0 +1,308 @@ +using ModelingToolkit, ModelingToolkitStandardLibrary, OrdinaryDiffEq +using ModelingToolkitStandardLibrary.Blocks + +@parameters t + +#= +Testing strategy: +The general strategy is to test systems using simple intputs where the solution +is known on closed form. For algebraic systems (without differential variables), +an integrator with a constant input is often used together with the system under test. +=# + +@testset "Constant" begin + @named c = Constant(; k=1) + @named int = Integrator() + @named iosys = ODESystem(connect(c.output, int.input), t, systems=[int, c]) + sys = structural_simplify(iosys) + prob = ODEProblem(sys, Pair[int.x=>1.0], (0.0, 1.0)) + sol = solve(prob, Rodas4()) + @test sol[int.output.u][end] ≈ 2 +end + +@testset "Derivative" begin + @named source = Sine(; frequency=1) + @named int = Integrator(; k=1) + @named der = Derivative(; k=1, T=0.001) + @named iosys = ODESystem([ + connect(source.output, der.input), + connect(der.output, int.input), + ], + t, + systems=[int, source, der], + ) + sys = structural_simplify(iosys) + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0)) + sol = solve(prob, Rodas4()) + @test isapprox(sol[source.output.u], sol[int.output.u], atol=1e-1) +end + +@testset "PT1" begin + pt1_func(t, k, T) = k * (1 - exp(-t / T)) # Known solution to first-order system + + k, T = 1.0, 0.1 + @named c = Constant(; k=1) + @named pt1 = FirstOrder(; k=k, T=T) + @named iosys = ODESystem(connect(c.output, pt1.input), t, systems=[pt1, c]) + sys = structural_simplify(iosys) + prob = ODEProblem(sys, Pair[], (0.0, 100.0)) + sol = solve(prob, Rodas4()) + @test sol[pt1.output.u] ≈ pt1_func.(sol.t, k, T) atol=1e-3 +end + +@testset "PT2" begin + # Known solution to second-order system + function pt2_func(t, k, w, d) + y = if d==0 + -k*(-1 + cos(t*w)) + else + d = complex(d) + real(k*(1 + (-cosh(sqrt(-1 + d^2)*t*w) - (d*sinh(sqrt(-1 + d^2)*t*w))/sqrt(-1 + d^2))/exp(d*t*w))) + end + end + + k, w, d = 1.0, 1.0, 0.5 + @named c = Constant(; k=1) + @named pt2 = SecondOrder(; k=k, w=w, d=d) + @named iosys = ODESystem(connect(c.output, pt2.input), t, systems=[pt2, c]) + sys = structural_simplify(iosys) + prob = ODEProblem(sys, Pair[], (0.0, 100.0)) + sol = solve(prob, Rodas4()) + @test sol[pt2.output.u] ≈ pt2_func.(sol.t, k, w, d) atol=1e-3 +end + +@testset "StateSpace" begin + A = [0 1;-1 -0.5] + B = [0, 1] + C = [0.9 1;] + D = [0;;] + @named ss = StateSpace(;A,B,C,D,x_start=zeros(2)) + @named c = Constant(; k=1) + @named model = ODESystem( + [ + connect(c.output, ss.input), + ], + t, + systems=[ss, c] + ) + sys = structural_simplify(model) + prob = ODEProblem(sys, Pair[], (0.0, 100.0)) + sol = solve(prob, Rodas4()) + # equilibrium point is at [1, 0] + @test sol[ss.x[1]][end] ≈ 1 atol=1e-3 + @test sol[ss.x[2]][end] ≈ 0 atol=1e-3 +end + +"""Second order demo plant""" +function Plant(;name, x_start=zeros(2)) + @named input = RealInput() + @named output = RealOutput() + D = Differential(t) + sts = @variables x1(t)=x_start[1] x2(t)=x_start[2] + eqs= [ + D(x1) ~ x2 + D(x2) ~ -x1 - 0.5 * x2 + input.u + output.u ~ 0.9 * x1 + x2 + ] + compose(ODESystem(eqs, t, sts, []; name), [input, output]) +end + +@testset "PI" begin + @named ref = Constant(; k=2) + @named pi_controller = PI(k=1, T=1) + @named plant = Plant() + @named fb = Feedback() + @named model = ODESystem( + [ + connect(ref.output, fb.input1), + connect(plant.output, fb.input2), + connect(fb.output, pi_controller.err_input), + connect(pi_controller.ctr_output, plant.input), + ], + t, + systems=[pi_controller, plant, ref, fb] + ) + sys = structural_simplify(model) + prob = ODEProblem(sys, Pair[], (0.0, 100.0)) + sol = solve(prob, Rodas4()) + @test sol[ref.output.u - plant.output.u][end] ≈ 0 atol=1e-3 # zero control error after 100s +end + +@testset "PID" begin + @named ref = Constant(; k=2) + @named pid_controller = PID(k=3, Ti=0.5, Td=100) + @named plant = Plant() + @named fb = Feedback() + @named model = ODESystem( + [ + connect(ref.output, fb.input1), + connect(plant.output, fb.input2), + connect(fb.output, pid_controller.err_input), + connect(pid_controller.ctr_output, plant.input), + ], + t, + systems=[pid_controller, plant, ref, fb] + ) + sys = structural_simplify(model) + prob = ODEProblem(sys, Pair[], (0.0, 100.0)) + sol = solve(prob, Rodas4()) + @test sol[ref.output.u - plant.output.u][end] ≈ 0 atol=1e-3 # zero control error after 100s +end + +@test_skip begin + k = 2 + Ti = 0.5 + Td = 0.7 + wp = 1 + wd = 1 + Ni = √(Td / Ti) + Nd = 12 + y_max = Inf + y_min = -Inf + u_r = sin(t) + u_y = 0 + function solve_with_input(; u_r, u_y, + controller = PID(; k, Ti, Td, wp, wd, Ni, Nd, y_max, y_min, name=:controller) + ) + @test count(ModelingToolkit.isinput, states(controller)) == 5 # 2 in PID, 1 sat, 1 I, 1 D + @test count(ModelingToolkit.isoutput, states(controller)) == 4 + # TODO: check number of unbound inputs when available, should be 2 + @named iosys = ODESystem([controller.u_r~u_r, controller.u_y~u_y], t, systems=[controller]) + sys = structural_simplify(iosys) + prob = ODEProblem(sys, Pair[], (0.0, 10.0)) + sol = solve(prob, Rodas4(), saveat=0:0.1:10) + controller, sys, sol + end + + # linearity in u_r + controller, sys, sol1 = solve_with_input(u_r=sin(t), u_y=0) + controller, sys, sol2 = solve_with_input(u_r=2sin(t), u_y=0) + @test sum(abs, sol1[controller.ea]) < eps() # This is the acutator model error due to saturation + @test 2sol1[controller.y] ≈ sol2[controller.y] rtol=1e-3 # linearity in u_r + + # linearity in u_y + controller, sys, sol1 = solve_with_input(u_y=sin(t), u_r=0) + controller, sys, sol2 = solve_with_input(u_y=2sin(t), u_r=0) + @test sum(abs, sol1[controller.ea]) < eps() # This is the acutator model error due to saturation + @test 2sol1[controller.y] ≈ sol2[controller.y] rtol=1e-3 # linearity in u_y + + # zero error + controller, sys, sol1 = solve_with_input(u_y=sin(t), u_r=sin(t)) + @test sum(abs, sol1[controller.y]) ≈ 0 atol=sqrt(eps()) + + # test saturation + controller, sys, sol1 = solve_with_input(; u_r=10sin(t), u_y=0, + controller = PID(; k, Ti, Td, wp, wd=0, Ni, Nd, y_max=10, y_min=-10, name=:controller) + ) + @test extrema(sol1[controller.y]) == (-10, 10) + + + # test P set-point weighting + controller, sys, sol1 = solve_with_input(; u_r=sin(t), u_y=0, + controller = PID(; k, Ti, Td, wp=0, wd, Ni, Nd, y_max, y_min, name=:controller) + ) + @test sum(abs, sol1[controller.ep]) ≈ 0 atol=sqrt(eps()) + + # test D set-point weighting + controller, sys, sol1 = solve_with_input(; u_r=sin(t), u_y=0, + controller = PID(; k, Ti, Td, wp, wd=0, Ni, Nd, y_max, y_min, name=:controller) + ) + @test sum(abs, sol1[controller.ed]) ≈ 0 atol=sqrt(eps()) + + + # zero integral gain + controller, sys, sol1 = solve_with_input(; u_r=sin(t), u_y=0, + controller = PID(; k, Ti=false, Td, wp, wd, Ni, Nd, y_max, y_min, name=:controller) + ) + @test isapprox(sum(abs, sol1[controller.I.y]), 0, atol=sqrt(eps())) + + + # zero derivative gain + @test_skip begin # During the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a non-finite number: [5] + controller, sys, sol1 = solve_with_input(; u_r=sin(t), u_y=0, + controller = PID(; k, Ti, Td=false, wp, wd, Ni, Nd, y_max, y_min, name=:controller) + ) + @test isapprox(sum(abs, sol1[controller.D.y]), 0, atol=sqrt(eps())) + end + + # Tests below can be activated when the concept of unbound_inputs exists in MTK + # @test isequal(Set(unbound_inputs(controller)), @nonamespace(Set([controller.u_r, controller.u_y]))) + # @test isempty(unbound_inputs(sys)) + # @test isequal(bound_inputs(sys), inputs(sys)) + # @test isequal( + # Set(bound_inputs(sys)), + # Set([controller.u_r, controller.u_y, controller.I.u, controller.D.u, controller.sat.u]) + # ) +end + +@testset "LimPI" begin + @named ref = Constant(; k=1) + @named pi_controller_lim = LimPI(k=3, T=0.5, u_max=1.5, u_min=-1.5, Ta=0.1) + @named pi_controller = PI(k=3, T=0.5) + @named sat = Limiter(y_max=1.5, y_min=-1.5) + @named plant = Plant() + @named fb = Feedback() + + # without anti-windup measure + sol = let + @named model = ODESystem( + [ + connect(ref.output, fb.input1), + connect(plant.output, fb.input2), + connect(fb.output, pi_controller.err_input), + connect(pi_controller.ctr_output, sat.input), + connect(sat.output, plant.input), + ], + t, + systems=[pi_controller, plant, ref, fb, sat] + ) + sys = structural_simplify(model) + prob = ODEProblem(sys, Pair[], (0.0, 20.0)) + sol = solve(prob, Rodas4()) + end + + # with anti-windup measure + sol_lim = let + @named model = ODESystem( + [ + connect(ref.output, fb.input1), + connect(plant.output, fb.input2), + connect(fb.output, pi_controller_lim.err_input), + connect(pi_controller_lim.ctr_output, sat.input), + connect(sat.output, plant.input), + ], + t, + systems=[pi_controller_lim, plant, ref, fb, sat] + ) + sys = structural_simplify(model) + prob = ODEProblem(sys, Pair[], (0.0, 20.0)) + sol = solve(prob, Rodas4()) + end + + @test sol[ref.output.u - plant.output.u][end] ≈ 0 atol=1e-3 # zero control error after 100s + @test sol_lim[ref.output.u - plant.output.u][end] ≈ 0 atol=1e-3 # zero control error after 100s + + # Plots.plot(sol; vars=[plant.output.u]) # without anti-windup measure + # Plots.plot!(sol_lim; vars=[plant.output.u]) # with anti-windup measure +end + +@testset "LimPID" begin + @named ref = Constant(; k=1) + @named pid_controller = LimPID(k=3, Ti=0.5, Td=100, u_max=1.5, u_min=-1.5, Ni=0.1/0.5) + @named plant = Plant() + @named model = ODESystem( + [ + connect(ref.output, pid_controller.reference), + connect(plant.output, pid_controller.measurement), + connect(pid_controller.ctr_output, plant.input), + ], + t, + systems=[pid_controller, plant, ref] + ) + sys = structural_simplify(model) + prob = ODEProblem(sys, Pair[], (0.0, 100.0)) + sol = solve(prob, Rodas4()) + + # Plots.plot(sol, vars=[plant.output.u, plant.input.u]) + @test sol[ref.output.u - plant.output.u][end] ≈ 0 atol=1e-3 # zero control error after 100s +end \ No newline at end of file diff --git a/test/Blocks/math.jl b/test/Blocks/math.jl new file mode 100644 index 000000000..c6bf82fe7 --- /dev/null +++ b/test/Blocks/math.jl @@ -0,0 +1,230 @@ +using ModelingToolkitStandardLibrary.Blocks +using ModelingToolkit, OrdinaryDiffEq + +@parameters t + + +@testset "Gain" begin + @named c = Constant(; k=1) + @named gain = Gain(1;) + @named int = Integrator(; k=1) + @named model = ODESystem([connect(c.output, gain.input), connect(gain.output, int.input)], t, systems=[int, gain, c]) + sys = structural_simplify(model) + + prob = ODEProblem(sys, Pair[int.x=>1.0], (0.0, 1.0)) + + sol = solve(prob, Rodas4()) + + @test sol[int.output.u][end] ≈ 2 +end + +@testset "Feedback loop" begin + @named c = Constant(; k=2) + @named gain = Gain(1;) + @named int = Integrator(; k=1) + @named fb = Feedback(;) + @named model = ODESystem( + [ + connect(c.output, fb.input1), + connect(fb.input2, int.output), + connect(fb.output, gain.input), + connect(gain.output, int.input), + ], + t, + systems=[int, gain, c, fb] + ) + sys = structural_simplify(model) + + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 100.0)) + + sol = solve(prob, Rodas4()) + @test sol[int.output.u][end] ≈ 2 +end + +@testset "Add" begin + @named c1 = Constant(; k=1) + @named c2 = Constant(; k=2) + @named add = Add(;) + @named int = Integrator(; k=1) + @named model = ODESystem( + [ + connect(c1.output, add.input1), + connect(c2.output, add.input2), + connect(add.output, int.input), + ], + t, + systems=[int, add, c1, c2] + ) + sys = structural_simplify(model) + + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0)) + + sol = solve(prob, Rodas4()) + @test sol[int.output.u][end] ≈ 3 +end + +@testset "Product" begin + @named c1 = Constant(; k=1) + @named c2 = Constant(; k=2) + @named prod = Product(;) + @named int = Integrator(; k=1) + @named model = ODESystem( + [ + connect(c1.output, prod.input1), + connect(c2.output, prod.input2), + connect(prod.output, int.input), + ], + t, + systems=[int, prod, c1, c2] + ) + sys = structural_simplify(model) + + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0)) + + sol = solve(prob, Rodas4()) + @test sol[int.output.u][end] ≈ 2 +end + +@testset "Division" begin + @named c1 = Constant(; k=1) + @named c2 = Constant(; k=2) + @named div = Division(;) + @named int = Integrator(; k=1) + @named model = ODESystem( + [ + connect(c1.output, div.input1), + connect(c2.output, div.input2), + connect(div.output, int.input), + ], + t, + systems=[int, div, c1, c2] + ) + sys = structural_simplify(model) + + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0)) + + sol = solve(prob, Rodas4()) + @test sol[int.output.u][end] ≈ 1/2 +end + +@testset "Abs" begin + @named c = Constant(; k=-1) + @named abs = Abs(;) + @named int = Integrator(; k=1) + @named model = ODESystem( + [ + connect(c.output, abs.input), + connect(abs.output, int.input), + ], + t, + systems=[int, abs, c] + ) + sys = structural_simplify(model) + + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0)) + + sol = solve(prob, Rodas4()) + @test sol[int.output.u][end] ≈ 1 +end + +@testset "Sqrt" begin + @named c = Constant(; k=4) + @named sqr = Sqrt(;) + @named int = Integrator(; k=1) + @named model = ODESystem( + [ + connect(c.output, sqr.input), + connect(sqr.output, int.input), + ], + t, + systems=[int, sqr, c] + ) + sys = structural_simplify(model) + + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0)) + + sol = solve(prob, Rodas4()) + @test sol[int.output.u][end] ≈ 2 +end + +@testset "Sign" begin + @named c = Constant(; k=3) + @named sig = Sign(;) + @named int = Integrator(; k=1) + @named model = ODESystem( + [ + connect(c.output, sig.input), + connect(sig.output, int.input), + ], + t, + systems=[int, sig, c] + ) + sys = structural_simplify(model) + + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0)) + + sol = solve(prob, Rodas4()) + @test sol[int.output.u][end] ≈ 1 +end + +@testset "MatrixGain" begin + K = [1 2; 3 4] + @named gain = MatrixGain(K;) + # TODO: +end + +@testset "Sum" begin + @named s = Sum(2;) + # TODO: +end + +@testset "Math" begin + for (block, func) in [(Abs, abs), (Sin, sin), (Cos, cos), (Tan, tan), (Asin, asin), (Acos, acos), (Atan, atan), (Sinh, sinh), (Cosh, cosh), (Tanh, tanh), (Exp, exp)] + @named source = Sine(frequency=1) + @named b = block() + @named int = Integrator() + @named model = ODESystem([connect(source.output, b.input), connect(b.output, int.input)], t, systems=[int, b, source]) + sys = structural_simplify(model) + + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0)) + + sol = solve(prob, Rodas4()) + @test sol[b.output.u] ≈ func.(sol[source.output.u]) + end + + # input must be positive + for (block, func) in [(Sqrt, sqrt), (Log, log), (Log10, log10)] + @named source = Sine(; frequency=1, offset=2) + @named b = block() + @named int = Integrator() + @named model = ODESystem([connect(source.output, b.input), connect(b.output, int.input)], t, systems=[int, b, source]) + sys = structural_simplify(model) + + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0)) + + sol = solve(prob, Rodas4()) + @test sol[b.output.u] ≈ func.(sol[source.output.u]) + end +end + +@testset "Atan2" begin + @named c1 = Constant(; k=1) + @named c2 = Constant(; k=2) + @named b = Atan2(;) + @named int = Integrator(; k=1) + @named model = ODESystem( + [ + connect(c1.output, b.input1), + connect(c2.output, b.input2), + connect(b.output, int.input), + ], + t, + systems=[int, b, c1, c2] + ) + sys = structural_simplify(model) + + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0)) + + sol = solve(prob, Rodas4()) + @test sol[int.output.u][end] ≈ atan(1, 2) +end diff --git a/test/Blocks/nonlinear.jl b/test/Blocks/nonlinear.jl new file mode 100644 index 000000000..2c792930e --- /dev/null +++ b/test/Blocks/nonlinear.jl @@ -0,0 +1,113 @@ +using ModelingToolkit, OrdinaryDiffEq +using ModelingToolkitStandardLibrary.Blocks +using ModelingToolkitStandardLibrary.Blocks: _clamp, _dead_zone + +@parameters t + +@testset "Limiter" begin + @testset "Constant" begin + @named c = Constant(; k=1) + @named int = Integrator(; k=1) + @named sat = Limiter(; y_min=-0.6, y_max=0.8) + @named model = ODESystem([ + connect(c.output, int.input), + connect(int.output, sat.input), + ], + t, + systems=[int, c, sat], + ) + sys = structural_simplify(model) + prob = ODEProblem(sys, [int.x=>1.0], (0.0, 1.0)) + + sol = solve(prob, Rodas4()) + @test sol[int.output.u][end] ≈ 2 + @test sol[sat.output.u][end] ≈ 0.8 + end + + @testset "Sine" begin + y_min, y_max = -0.3, 0.5 + @named source = Sine(; frequency=1/2) + @named lim = Limiter(; y_max=y_max, y_min=y_min) + @named int = Integrator(; k=1) + @named iosys = ODESystem([ + connect(source.output, lim.input), + connect(lim.output, int.input), + ], + t, + systems=[source, lim, int], + ) + sys = structural_simplify(iosys) + + prob = ODEProblem(sys, Pair[], (0.0, 10.0)) + + sol = solve(prob, Rodas4()) + @test all(abs.(sol[lim.output.u]) .<= 0.5) + @test all(isapprox.(sol[lim.output.u], _clamp.(sol[source.output.u], y_min, y_max), atol=1e-2)) + + # Plots.plot(sol; vars=[source.output.u, lim.output.u]) + # Plots.scatter(sol[source.output.u], sol[lim.output.u]) + # Plots.scatter!(sol[source.output.u], _clamp.(sol[source.output.u], y_min, y_max)) + end +end + +@testset "DeadZone" begin + @testset "Constant" begin + @named c = Constant(; k=1) + @named int = Integrator(; k=1) + @named dz = DeadZone(; u_min=-2, u_max=1) + @named model = ODESystem([ + connect(c.output, int.input), + connect(int.output, dz.input), + ], + t, + systems=[int, c, dz], + ) + sys = structural_simplify(model) + prob = ODEProblem(sys, [int.x=>1.0], (0.0, 1.0)) + sol = solve(prob, Rodas4()) + + @test sol[int.output.u][end] ≈ 2 + end + + @testset "Sine" begin + u_min, u_max = -2, 1 + @named source = Sine(; amplitude=3, frequency=1/2) + @named dz = DeadZone(; u_min=u_min, u_max=u_max) + @named int = Integrator(; k=1) + @named model = ODESystem([ + connect(source.output, dz.input), + connect(dz.output, int.input), + ], + t, + systems=[int, source, dz], + ) + sys = structural_simplify(model) + prob = ODEProblem(sys, [int.x=>1.0], (0.0, 10.0)) + sol = solve(prob, Rodas4()) + + @test all(sol[dz.output.u] .<= 2) + @test all(sol[dz.output.u] .>= -1) + @test all(isapprox.(sol[dz.output.u], _dead_zone.(sol[source.output.u], u_min, u_max), atol=1e-2)) + + # Plots.plot(sol; vars=[source.output.u, dz.output.u]) + # Plots.scatter(sol[source.output.u], sol[dz.output.u]) + # Plots.scatter!(sol[source.output.u], _dead_zone.(sol[source.output.u], u_min, u_max)) + end +end + +@testset "SlewRateLimiter" begin + @named source = Sine(; frequency=1/2) + @named rl = SlewRateLimiter(; rising=1, falling=-1, Td=0.001, y_start=-1/3) + @named iosys = ODESystem([ + connect(source.output, rl.input), + ], + t, + systems=[source, rl], + ) + sys = structural_simplify(iosys) + + prob = ODEProblem(sys, Pair[], (0.0, 10.0)) + + sol = solve(prob, Rodas4(), saveat=0.01, abstol=1e-10, reltol=1e-10) + @test all(abs.(sol[rl.output.u]) .<= 0.51) +end \ No newline at end of file diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl new file mode 100644 index 000000000..2ac7f22cf --- /dev/null +++ b/test/Blocks/sources.jl @@ -0,0 +1,155 @@ +using ModelingToolkit, ModelingToolkitStandardLibrary, OrdinaryDiffEq +using ModelingToolkitStandardLibrary.Blocks + +@parameters t + +@testset "Constant" begin + @named src = Constant(k=2) + @named int = Integrator() + @named iosys = ODESystem([ + connect(src.output, int.input), + ], + t, + systems=[int, src], + ) + sys = structural_simplify(iosys) + + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0)) + + sol = solve(prob, Rodas4()) + @test sol[src.output.u][end] ≈ 2 atol=1e-3 +end + +@testset "Sine" begin + sine(t, frequency, amplitude, phase, offset, start_time) = offset + ifelse(t < start_time, 0, amplitude* sin(2*pi*frequency*(t - start_time) + phase)) + + frequency=1 + amplitude=2 + phase=0 + offset=1 + start_time=0 + + @named src = Sine(frequency=frequency, amplitude=amplitude, phase=phase, offset=offset, start_time=start_time) + @named int = Integrator() + @named iosys = ODESystem([ + connect(src.output, int.input), + ], + t, + systems=[int, src], + ) + sys = structural_simplify(iosys) + + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0)) + + sol = solve(prob, Rodas4()) + @test sol[src.output.u] ≈ sine.(sol.t, frequency, amplitude, phase, offset, start_time) atol=1e-3 +end + +@testset "Cosine" begin + cosine(t, frequency, amplitude, phase, offset, start_time) = offset + ifelse(t < start_time, 0, amplitude* cos(2*pi*frequency*(t - start_time) + phase)) + + frequency=1 + amplitude=2 + phase=0 + offset=1 + start_time=0 + + @named src = Cosine(frequency=frequency, amplitude=amplitude, phase=phase, offset=offset, start_time=start_time) + @named int = Integrator() + @named iosys = ODESystem([ + connect(src.output, int.input), + ], + t, + systems=[int, src], + ) + sys = structural_simplify(iosys) + + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0)) + + sol = solve(prob, Rodas4()) + @test sol[src.output.u] ≈ cosine.(sol.t, frequency, amplitude, phase, offset, start_time) atol=1e-3 +end + +@testset "ContinuousClock" begin + cont_clock(t, offset, start_time) = offset + ifelse(t < start_time, 0, t - start_time) + + offset, start_time = 1, 0 + + @named src = ContinuousClock(offset=offset, start_time=start_time) + @named int = Integrator() + @named iosys = ODESystem([ + connect(src.output, int.input), + ], + t, + systems=[int, src], + ) + sys = structural_simplify(iosys) + + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0)) + + sol = solve(prob, Rodas4()) + @test sol[src.output.u] ≈ cont_clock.(sol.t, offset, start_time) atol=1e-3 +end + +@testset "Ramp" begin + ramp(t, offset, height, duration, start_time) = offset + ifelse(t < start_time, 0, ifelse(t < (start_time + duration), (t - start_time) * height / duration, height)) + + offset, height, duration, start_time = 1, 2, 2, 0 + + @named src = Ramp(offset=offset, height=height, duration=duration, start_time=start_time) + @named int = Integrator() + @named iosys = ODESystem([ + connect(src.output, int.input), + ], + t, + systems=[int, src], + ) + sys = structural_simplify(iosys) + + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0)) + + sol = solve(prob, Rodas4()) + @test sol[src.output.u] ≈ ramp.(sol.t, offset, height, duration, start_time) atol=1e-3 +end + +@testset "Step" begin + step(t, offset, height, start_time) = offset + ifelse(t < start_time, 0, height) + + offset, height, start_time = 1, 2, 5 + + @named src = Step(offset=offset, height=height, start_time=start_time) + @named int = Integrator() + @named iosys = ODESystem([ + connect(src.output, int.input), + ], + t, + systems=[int, src], + ) + sys = structural_simplify(iosys) + + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0)) + + sol = solve(prob, Rodas4()) + @test sol[src.output.u] ≈ step.(sol.t, offset, height, start_time) atol=1e-3 +end + +@testset "ExpSine" begin + exp_sine(t, amplitude, frequency, damping, phase, start_time) = offset + ifelse(t < start_time, 0, amplitude * exp(-damping * (t - start_time)) * sin(2*pi*frequency*(t - start_time) + phase)) + + frequency, amplitude, damping, phase, offset, start_time = 3, 2, 0.10, 0, 0, 0 + + @named src = ExpSine(frequency=frequency, amplitude=amplitude, damping=damping, phase=phase, offset=offset, start_time=start_time) + @named int = Integrator() + @named iosys = ODESystem([ + connect(src.output, int.input), + ], + t, + systems=[int, src], + ) + sys = structural_simplify(iosys) + + prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0)) + + sol = solve(prob, Rodas4()) + @test sol[src.output.u] ≈ exp_sine.(sol.t, amplitude, frequency, damping, phase, start_time) atol=1e-3 +end \ No newline at end of file diff --git a/test/Electrical/analog.jl b/test/Electrical/analog.jl new file mode 100644 index 000000000..fc7184d72 --- /dev/null +++ b/test/Electrical/analog.jl @@ -0,0 +1,333 @@ +using ModelingToolkitStandardLibrary.Electrical, ModelingToolkit, OrdinaryDiffEq, Test +using ModelingToolkitStandardLibrary.Electrical: _step, _square_wave, _triangular_wave, _cos_wave, _damped_sine_wave, _ramp + +# using Plots + +@parameters t + +@testset "sensors" begin + @named source = ConstantVoltage(V=10) + @named resistor = Resistor(R=1) + @named capacitor = Capacitor(C=1) + @named ground = Ground() + + @named voltage_sensor = VoltageSensor() + @named current_sensor = CurrentSensor() + @named power_sensor = PowerSensor() + + connections = [ + connect(source.p, resistor.p) + connect(resistor.n, current_sensor.p) + connect(current_sensor.n, power_sensor.pc) + connect(power_sensor.nc, capacitor.p) + connect(capacitor.n, source.n, ground.g) + connect(capacitor.p, voltage_sensor.p) + connect(capacitor.n, voltage_sensor.n) + connect(capacitor.p, power_sensor.pv) + connect(capacitor.n, power_sensor.nv) + ] + + @named model = ODESystem(connections, t; systems=[resistor, capacitor, source, ground, voltage_sensor, current_sensor, power_sensor]) + sys = structural_simplify(model) + prob = ODAEProblem(sys, Pair[], (0.0, 10.0)) + sol = solve(prob, Tsit5()) + + # Plots.plot(sol; vars=[capacitor.v, voltage_sensor.v]) + # Plots.plot(sol; vars=[power_sensor.power, capacitor.i * capacitor.v]) + # Plots.plot(sol; vars=[resistor.i, current_sensor.i]) + @test sol[capacitor.v] ≈ sol[voltage_sensor.v] atol=1e-3 + @test sol[power_sensor.power] ≈ sol[capacitor.i * capacitor.v] atol=1e-3 + @test sol[resistor.i] ≈ sol[current_sensor.i] atol=1e-3 +end + +# simple voltage divider +@testset "voltage divider" begin + @named source = ConstantVoltage(V=10) + @named R1 = Resistor(R=1e3) + @named R2 = Resistor(R=1e3) + @named ground = Ground() + + connections = [ + connect(source.p, R1.p) + connect(R1.n, R2.p) + connect(R2.n, source.n, ground.g) + ] + + @named model = ODESystem(connections, t, systems=[R1, R2, source, ground]) + sys = structural_simplify(model) + prob = ODEProblem(sys, Pair[], (0, 2.0)) + sol = solve(prob, Rodas4()) # has no state; does not work with Tsit5 + @test sol[R1.p.v][end] ≈ 10 atol=1e-3 + @test sol[R1.n.v][end] ≈ 5 atol=1e-3 + @test sol[R2.n.v][end] ≈ 0 atol=1e-3 +end + +# simple RC +@testset "RC" begin + @named source = ConstantVoltage(V=10) + @named resistor = Resistor(R=1) + @named capacitor = Capacitor(C=1) + @named ground = Ground() + + connections = [ + connect(source.p, resistor.p) + connect(resistor.n, capacitor.p) + connect(capacitor.n, source.n, ground.g) + ] + + @named model = ODESystem(connections, t; systems=[resistor, capacitor, source, ground]) + sys = structural_simplify(model) + prob = ODAEProblem(sys, [capacitor.v => 0.0], (0.0, 10.0)) + sol = solve(prob, Tsit5()) + + # Plots.plot(sol; vars=[source.v, capacitor.v]) + @test sol[capacitor.v][end] ≈ 10 atol=1e-3 +end + +# simple RL +@testset "RL" begin + @named source = ConstantVoltage(V=10) + @named resistor = Resistor(R=1) + @named inductor = Inductor(L=1.0) + @named ground = Ground() + + connections = [ + connect(source.p, resistor.p) + connect(resistor.n, inductor.p) + connect(inductor.n, source.n, ground.g) + ] + + @named model = ODESystem(connections, t; systems=[resistor, inductor, source, ground]) + sys = structural_simplify(model) + prob = ODAEProblem(sys, [inductor.i => 0.0], (0.0, 10.0)) + sol = solve(prob, Tsit5()) + + # Plots.plot(sol; vars=[inductor.i, inductor.i]) + @test sol[inductor.i][end] ≈ 10 atol=1e-3 +end + +# RC with different voltage sources +@testset "RC with voltage sources" begin + @named source_const = ConstantVoltage(V=10) + @named source_sin = SineVoltage(offset=1, amplitude=10, frequency=2, start_time=0.5, phase=0) + @named source_step = StepVoltage(offset=1, height=10, start_time=0.5) + @named source_tri = TriangularVoltage(offset=1, start_time=0.5, amplitude=10, frequency=2) + @named source_dsin = ExpSineVoltage(offset=1, amplitude=10, frequency=2, start_time=0.5, phase=0, damping=0.5) + @named source_ramp = RampVoltage(offset=1, height=10, start_time=0.5, duration=1) + sources = [source_const, source_sin, source_step, source_tri, source_dsin, source_ramp] + + @named resistor = Resistor(R=1) + @named capacitor = Capacitor(C=1) + @named ground = Ground() + + for source in sources + connections = [ + connect(source.p, resistor.p) + connect(resistor.n, capacitor.p) + connect(capacitor.n, source.n, ground.g) + ] + + @named model = ODESystem(connections, t; systems=[resistor, capacitor, source, ground]) + sys = structural_simplify(model) + prob = ODAEProblem(sys, [capacitor.v => 0.0], (0.0, 10.0)) + @test_nowarn sol = solve(prob, Tsit5()) + @test_nowarn sol = solve(prob, Rodas4()) + + # Plots.plot(sol; vars=[source.v, capacitor.v]) + end +end + +# RL with different voltage sources +@testset "RL with voltage sources" begin + @named source_const = ConstantVoltage(V=10) + @named source_sin = SineVoltage(offset=1, amplitude=10, frequency=2, start_time=0.5, phase=0) + @named source_step = StepVoltage(offset=1, height=10, start_time=0.5) + @named source_tri = TriangularVoltage(offset=1, start_time=0.5, amplitude=10, frequency=2) + @named source_dsin = ExpSineVoltage(offset=1, amplitude=10, frequency=2, start_time=0.5, phase=0, damping=0.5) + @named source_ramp = RampVoltage(offset=1, height=10, start_time=0.5, duration=1) + sources = [source_const, source_sin, source_step, source_tri, source_dsin, source_ramp] + + @named resistor = Resistor(R=1.0) + @named inductor = Inductor(L=1.0) + @named ground = Ground() + + for source in sources + connections = [ + connect(source.p, resistor.p) + connect(resistor.n, inductor.p) + connect(inductor.n, source.n, ground.g) + ] + + @named model = ODESystem(connections, t; systems=[resistor, inductor, source, ground]) + sys = structural_simplify(model) + prob = ODAEProblem(sys, [inductor.i => 0.0], (0.0, 10.0)) + @test_nowarn sol = solve(prob, Tsit5()) + @test_nowarn sol = solve(prob, Rodas4()) + + # Plots.plot(sol; vars=[source.i, inductor.i]) + end +end + +# RC with different current sources +@testset "RC with current sources" begin + @named source_const = ConstantCurrent(I=10) + @named source_sin = SineCurrent(offset=1, amplitude=10, frequency=2, start_time=0.5, phase=0) + @named source_step = StepCurrent(offset=1, height=10, start_time=0.5) + @named source_tri = TriangularCurrent(offset=1, start_time=0.5, amplitude=10, frequency=2) + @named source_dsin = ExpSineCurrent(offset=1, amplitude=10, frequency=2, start_time=0.5, phase=0, damping=0.5) + @named source_ramp = RampCurrent(offset=1, height=10, start_time=0.5, duration=1) + sources = [source_const, source_sin, source_step, source_tri, source_dsin, source_ramp] + + @named resistor = Resistor(R=1) + @named capacitor = Capacitor(C=1) + @named ground = Ground() + + for source in sources + connections = [ + connect(source.p, resistor.p) + connect(resistor.n, capacitor.p) + connect(capacitor.n, source.n, ground.g) + ] + + @named model = ODESystem(connections, t; systems=[resistor, capacitor, source, ground]) + sys = structural_simplify(model) + prob = ODAEProblem(sys, [capacitor.v => 0.0], (0.0, 10.0)) + @test_nowarn sol = solve(prob, Tsit5()) + @test_broken sol = solve(prob, Rodas4()) + + # Plots.plot(sol; vars=[source.v, capacitor.v]) + end +end + +@testset "Integrator" begin + R=1e3 + f=1 + Vin=5 + @named ground = Ground() + @named R1 = Resistor(R=R) + @named R2 = Resistor(R=100*R) + @named C1 = Capacitor(C=1/(2 * pi * f * R)) + @named opamp = IdealOpAmp() + @named square = SquareVoltage(amplitude=Vin) + @named sensor = VoltageSensor() + + connections = [ + connect(square.p, R1.p) + connect(R1.n, C1.n, R2.p, opamp.n1) + connect(opamp.p2, C1.p, R2.n) + connect(opamp.p1, ground.g, opamp.n2, square.n) + connect(opamp.p2, sensor.p) + connect(sensor.n, ground.g) + ] + @named model = ODESystem(connections, t, systems = [R1, R2, opamp, square, C1, ground, sensor]) + sys = structural_simplify(model) + u0 = [ + C1.v => 0.0 + R1.v => 0.0 + ] + prob = ODEProblem(sys, u0, (0, 100.0)) + sol = solve(prob, Rodas4()) + @test sol[opamp.v2] == sol[C1.v] # Not a great one however. Rely on the plot + @test sol[opamp.p2.v] == sol[sensor.v] + + # plot(sol, vars=[sensor.v, square.v, C1.v]) +end + +@testset "Voltage function generators" begin + st, o, h, f, A, et, ϕ, d, δ = 0.7, 1.25, 3, 2, 2.5, 2.5, π/4, 0.1, 0.0001 + + @named res = Resistor(R=1) + @named cap = Capacitor(C=1) + @named ground = Ground() + @named voltage_sensor = VoltageSensor() + @named vstep = StepVoltage(start_time=st, offset=o, height=h) + @named vsquare = SquareVoltage(offset=o, start_time=st, amplitude=A, frequency=f) + @named vtri = TriangularVoltage(offset=o, start_time=st, amplitude=A, frequency=f) + # @named vsawtooth = SawToothVoltage(amplitude=A, start_time=st, frequency=f, offset=o) + @named vcosine = CosineVoltage(offset=o, amplitude=A, frequency=f, start_time=st, phase=ϕ) + @named vdamped_sine = ExpSineVoltage(offset=o, amplitude=A, frequency=f, start_time=st, phase=ϕ, damping=d) + @named vramp = RampVoltage(offset=o, start_time=st, duration=et-st, height=h) + + vsources = [vtri, vsquare, vstep, vcosine, vdamped_sine, vramp] + waveforms(i, x) = getindex([o .+ (x .> st) .* _triangular_wave.(x, δ, f, A, st), + o .+ (x .> st) .* _square_wave.(x, δ, f, A, st), + o .+ _step.(x, δ, h, st), + # o .+ (x .> st). * _sawtooth_wave.(x, δ, f, A, st), + o .+ (x .> st) .* _cos_wave.(x, f, A, st, ϕ), + o .+ (x .> st) .* _damped_sine_wave.(x, f, A, st, ϕ, d), + o .+ _ramp.(x, δ, st, et, h)], i) + for i in 1:length(vsources) + vsource = vsources[i] + # @info Symbolics.getname(vsource) + eqs = [ + connect(vsource.p, voltage_sensor.p, res.p) + connect(res.n, cap.p) + connect(ground.g, voltage_sensor.n, vsource.n, cap.n) + ] + @named vmodel = ODESystem(eqs, t, systems = [voltage_sensor, res, cap, vsource, ground]) + vsys = structural_simplify(vmodel) + + u0 = [ + vsource.v => 1 + res.v => 1 + ] + + prob = ODAEProblem(vsys, u0, (0, 10.0)) + sol = solve(prob, dt=0.1, Tsit5()) + + @test sol[vsource.v][1150:end] ≈ waveforms(i, sol.t)[1150:end] atol=1e-1 + # For visual inspection + # plt = plot(sol; vars=[vsource.v]) + # savefig(plt, "test_voltage_$(Symbolics.getname(vsource))") + end +end + +@testset "Current function generators" begin + st, o, h, f, A, et, ϕ, d, δ = 0.7, 1.25, 3, 2, 2.5, 2.5, π/4, 0.1, 0.0001 + + @named ground = Ground() + @named res = Resistor(R=1.0) + @named cap = Capacitor(C=1) + @named current_sensor = CurrentSensor() + @named istep = StepCurrent(start_time=st, offset=o, height=h) + @named isquare = SquareCurrent(offset=o, start_time=st, amplitude=A, frequency=f) + @named itri = TriangularCurrent(offset=o, start_time=st, amplitude=A, frequency=f) + # @named isawtooth = SawToothCurrent(amplitude=A, start_time=st, frequency=f, offset=o) + @named icosine = CosineCurrent(offset=o, amplitude=A, frequency=f, start_time=st, phase=ϕ) + @named idamped_sine = ExpSineCurrent(offset=o, amplitude=A, frequency=f, start_time=st, phase=ϕ, damping=d) + @named iramp = RampCurrent(offset=o, start_time=st, duration=et-st, height=h) + + isources = [itri, isquare, istep, icosine, idamped_sine, iramp] + waveforms(i, x) = getindex([o .+ (x .> st) .* _triangular_wave.(x, δ, f, A, st), + o .+ (x .> st) .* _square_wave.(x, δ, f, A, st), + o .+ _step.(x, δ, h, st), + # o .+ (x .> st). * _sawtooth_wave.(x, δ, f, A, st), + o .+ (x .> st) .* _cos_wave.(x, f, A, st, ϕ), + o .+ (x .> st) .* _damped_sine_wave.(x, f, A, st, ϕ, d), + o .+ _ramp.(x, δ, st, et, h)], i) + + for i in 1:length(isources) + isource = isources[i] + eqs = [ + connect(isource.p, current_sensor.n) + connect(current_sensor.p, res.p) + connect(res.n, cap.p) + connect(isource.n, ground.g, cap.n) + ] + @named model = ODESystem(eqs, t, systems = [current_sensor, isource, res, cap, ground]) + isys = structural_simplify(model) + + u0 = [ + isource.i => 1.0 + res.v => 1.0 + cap.v => 0.0 + ] + prob = ODAEProblem(isys, u0, (0, 10.0)) + sol = solve(prob, Tsit5()) + + @test sol[isource.i][1150:end] ≈ waveforms(i, sol.t)[1150:end] atol=1e-1 + # For visual inspection + # plt = plot(sol) + # savefig(plt, "test_current_$(Symbolics.getname(isource))") + end +end \ No newline at end of file diff --git a/test/digital.jl b/test/Electrical/digital.jl similarity index 100% rename from test/digital.jl rename to test/Electrical/digital.jl diff --git a/test/Magnetic/magnetic.jl b/test/Magnetic/magnetic.jl new file mode 100644 index 000000000..7de2ad712 --- /dev/null +++ b/test/Magnetic/magnetic.jl @@ -0,0 +1,48 @@ +using ModelingToolkitStandardLibrary.Magnetic, ModelingToolkit, OrdinaryDiffEq, Test + +import ModelingToolkitStandardLibrary.Electrical +import ModelingToolkitStandardLibrary.Magnetic +using ModelingToolkit, OrdinaryDiffEq, Test +# using Plots + +@parameters t + +@testset "Inductor" begin + mu_air = 1 + l_air = 0.0001 + mu_Fe = 1000 + l_Fe = 4*0.065 + a = b = 0.25 + + @named source = Electrical.SineVoltage(amplitude=230*sqrt(2), frequency=50, phase=pi/2) + @named r = Electrical.Resistor(R=7.5) + @named ground = Electrical.Ground() + @named coil = Magnetic.FluxTubes.ElectroMagneticConverter(N=600) + @named ground_m = Magnetic.FluxTubes.Ground() + @named r_mAirPar = Magnetic.FluxTubes.ConstantReluctance(R_m=a * b * l_air * mu_air) + @named r_mFe = Magnetic.FluxTubes.ConstantReluctance(R_m=a * b * l_Fe * mu_Fe) + @named r_mLeak = Magnetic.FluxTubes.ConstantReluctance(R_m=1.2e6) + connections = [ + connect(source.p, r.p) + connect(r.n, coil.p) + connect(source.n, coil.n) + connect(coil.port_p, r_mLeak.port_p) + connect(r_mLeak.port_p, r_mAirPar.port_p) + connect(r_mAirPar.port_n, r_mFe.port_p) + connect(r_mFe.port_n, r_mLeak.port_n) + connect(r_mFe.port_n, coil.port_n) + connect(ground.g, source.n) + connect(ground_m.port, r_mFe.port_n) + ] + @named model = ODESystem(connections, t, systems=[source, r, ground, coil, ground_m, r_mAirPar, r_mFe, r_mLeak]) + sys = structural_simplify(model) + prob = ODEProblem(sys, Pair[], (0, 0.1)) + sol = solve(prob, Rodas4()) + + # Plots.plot(sol; vars=[r.i]) + # Plots.plot(sol; vars=[r_mFe.V_m, r_mFe.Phi]) + + @test sol[r_mFe.Phi] == sol[r_mAirPar.Phi] + @test all(sol[coil.port_p.Phi] + sol[r_mLeak.Phi] + sol[r_mAirPar.Phi] .== 0) +end + diff --git a/test/Mechanical/rotational.jl b/test/Mechanical/rotational.jl new file mode 100644 index 000000000..f1bbeb07d --- /dev/null +++ b/test/Mechanical/rotational.jl @@ -0,0 +1,101 @@ +using ModelingToolkitStandardLibrary.Mechanical.Rotational, ModelingToolkit, OrdinaryDiffEq, Test +import ModelingToolkitStandardLibrary.Blocks +# using Plots + +@parameters t + +@testset "two inertias" begin + @named fixed = Fixed() + @named inertia1 = Inertia(J=2) # this one is fixed + @named spring = Spring(c=1e4) + @named damper = Damper(d=10) + @named inertia2 = Inertia(J=2, phi_start=pi/2) + + connections = [ + connect(fixed.flange, inertia1.flange_b) + connect(inertia1.flange_b, spring.flange_a, damper.flange_a) + connect(spring.flange_b, damper.flange_b, inertia2.flange_a) + ] + + @named model = ODESystem(connections, t, systems=[fixed, inertia1, inertia2, spring, damper]) + sys = structural_simplify(model) + prob = ODEProblem(sys, Pair[], (0, 10.0)) + sol = solve(prob, Rodas4()) + + # Plots.plot(sol; vars=[inertia1.w, inertia2.w]) + + @test all(sol[inertia1.w] .== 0) + @test sol[inertia2.w][end] ≈ 0 atol=1e-3 # all energy has dissipated +end + +@testset "two inertias with driving torque" begin + amplitude=10 # Amplitude of driving torque + frequency=5 # Frequency of driving torque + J_motor=0.1 # Motor inertia + + @named fixed = Fixed() + @named torque = Torque(use_support=true) + @named inertia1 = Inertia(J=2, phi_start=pi/2) + @named spring = Spring(c=1e4) + @named damper = Damper(d=10) + @named inertia2 = Inertia(J=4) + @named sine = Blocks.Sine(amplitude=amplitude, frequency=frequency) + + connections = [ + connect(sine.output, torque.tau) + connect(torque.support, fixed.flange) + connect(torque.flange, inertia1.flange_a) + connect(inertia1.flange_b, spring.flange_a, damper.flange_a) + connect(spring.flange_b, damper.flange_b, inertia2.flange_a) + ] + + @named model = ODESystem(connections, t, systems=[fixed, torque, inertia1, inertia2, spring, damper, sine]) + sys = structural_simplify(model) + prob = ODAEProblem(sys, Pair[], (0, 1.0)) + sol = solve(prob, Rodas4()) + + # Plots.plot(sol; vars=[inertia1.w, -inertia2.w*2]) + + @test all(isapprox.(sol[inertia1.w], -sol[inertia2.w]*2, atol=1)) # exact opposite oscillation with smaller amplitude J2 = 2*J1 + @test all(sol[torque.flange.tau] .== -sol[sine.output.u]) # torque source is equal to negative sine +end + +# see: https://doc.modelica.org/Modelica%204.0.0/Resources/helpWSM/Modelica/Modelica.Mechanics.Rotational.Examples.First.html +@testset "first example" begin + amplitude=10 # Amplitude of driving torque + frequency=5 # Frequency of driving torque + J_motor=0.1 # Motor inertia + J_load=2 # Load inertia + ratio=10 # Gear ratio + damping=10 # Damping in bearing of gear + + @named fixed = Fixed() + @named torque = Torque(use_support=true) + @named inertia1 = Inertia(J=J_motor) + @named idealGear = IdealGear(ratio=ratio, use_support=true) + @named inertia2 = Inertia(J=2) + @named spring = Spring(c=1e4) + @named inertia3 = Inertia(J=J_load) + @named damper = Damper(d=damping) + @named sine = Blocks.Sine(amplitude=amplitude, frequency=frequency) + + connections = [ + connect(inertia1.flange_b, idealGear.flange_a) + connect(idealGear.flange_b, inertia2.flange_a) + connect(inertia2.flange_b, spring.flange_a) + connect(spring.flange_b, inertia3.flange_a) + connect(damper.flange_a, inertia2.flange_b) + connect(damper.flange_b, fixed.flange) + connect(sine.output, torque.tau) + connect(torque.support, fixed.flange) + connect(idealGear.support, fixed.flange) + connect(torque.flange, inertia1.flange_a) + ] + + @named model = ODESystem(connections, t, systems=[fixed, torque, inertia1, idealGear, inertia2, spring, inertia3, damper, sine]) + sys = structural_simplify(model) + @test_broken prob = ODAEProblem(sys, Pair[], (0, 1.0)) # KeyError: key 25 not found + # sol = solve(prob, Rodas4()) + + # Plots.plot(sol; vars=[inertia2.w, inertia3.w]) +end \ No newline at end of file diff --git a/test/Thermal/demo.jl b/test/Thermal/demo.jl new file mode 100644 index 000000000..8db93df4d --- /dev/null +++ b/test/Thermal/demo.jl @@ -0,0 +1,23 @@ +using ModelingToolkitStandardLibrary.Thermal, ModelingToolkit, OrdinaryDiffEq, Test +@parameters t + +# Modelica example +begin + @named mass1 = HeatCapacitor(C=15, T_start=373.15) + @named mass2 = HeatCapacitor(C=15, T_start=273.15) + @named conduction = ThermalConductor(G=10) + @named Tsensor1 = TemperatureSensor() + @named Tsensor2 = TemperatureSensor() + + connections = [ + connect(mass1.port, conduction.port_a), + connect(conduction.port_b, mass2.port), + connect(mass1.port, Tsensor1.port), + connect(mass2.port, Tsensor2.port), + ] + + @named model = ODESystem(connections, t, systems=[mass1, mass2, conduction, Tsensor1, Tsensor2]) + sys = structural_simplify(model) + prob = ODEProblem(sys, Pair[], (0, 3.0)) + sol = solve(prob, Rodas4()) +end \ No newline at end of file diff --git a/test/thermal.jl b/test/Thermal/thermal.jl similarity index 65% rename from test/thermal.jl rename to test/Thermal/thermal.jl index 47509492f..be95a7e16 100644 --- a/test/thermal.jl +++ b/test/Thermal/thermal.jl @@ -15,29 +15,27 @@ using ModelingToolkitStandardLibrary.Thermal, ModelingToolkit, OrdinaryDiffEq, T @info "Building a single-body system..." eqs = [ - connect(mass1.hp, th_conductor.hp1) - connect(th_conductor.hp2, reltem_sensor.hp1) - connect(reltem_sensor.hp2, tem_src.hp) + connect(mass1.port, th_conductor.port_a) + connect(th_conductor.port_b, reltem_sensor.port_a) + connect(reltem_sensor.port_b, tem_src.port) ] @named h1 = ODESystem(eqs, t, systems=[mass1, reltem_sensor, tem_src, th_conductor]) sys = structural_simplify(h1) u0 = [ - mass1.T => 2.0 - th_conductor.T => 10.0 + mass1.T => 2.0 ] prob = ODEProblem(sys, u0, (0, 2.0)) - sol = solve(prob, Rosenbrock23()) + sol = solve(prob, Rodas4()) - temperatures = reduce(hcat, sol.u) # Check if Relative temperature sensor reads the temperature of heat capacitor # when connected to a thermal conductor and a fixed temperature source - @test sol[reltem_sensor.T] == temperatures[1, :] - temperatures[2, :] - sol[tem_src.hp.T] + @test sol[reltem_sensor.T] + sol[tem_src.port.T] == sol[mass1.T] + sol[th_conductor.dT] @info "Building a two-body system..." eqs = [ - connect(T_sensor1.hp, mass1.hp, th_conductor.hp1) - connect(th_conductor.hp2, mass2.hp, T_sensor2.hp) + connect(T_sensor1.port, mass1.port, th_conductor.port_a) + connect(th_conductor.port_b, mass2.port, T_sensor2.port) final_T ~ (mass1.C * mass1.T + mass2.C * mass2.T) / (mass1.C + mass2.C) ] @@ -51,7 +49,7 @@ using ModelingToolkitStandardLibrary.Thermal, ModelingToolkit, OrdinaryDiffEq, T final_T => 12 ] prob = ODEProblem(sys, u0, (0, 3.0)) - sol = solve(prob, Rosenbrock23()) + sol = solve(prob, Rodas4()) m1, m2 = sol.u[end] @test m1 ≈ m2 atol=1e-1 @@ -60,22 +58,22 @@ using ModelingToolkitStandardLibrary.Thermal, ModelingToolkit, OrdinaryDiffEq, T @test sol[T_sensor2.T] == mass_T[2, :] end -# Test HeatFlowSensor, FixedHeatFlow, ThermalResistor, ThermalConductor, ThermalGround +# Test HeatFlowSensor, FixedHeatFlow, ThermalResistor, ThermalConductor @testset "Heat flow system" begin C, G, R = 10, 10, 10 - @named flow_src = FixedHeatFlow(Q_flow=50, α=100) + @named flow_src = FixedHeatFlow(Q_flow=50, alpha=100) @named mass1 = HeatCapacitor(C=C) @named hf_sensor1 = HeatFlowSensor() @named hf_sensor2 = HeatFlowSensor() @named th_conductor = ThermalConductor(G=G) @named th_resistor = ThermalResistor(R=R) - @named th_ground = ThermalGround() + @named th_ground = FixedTemperature(T=0) @info "Building a heat-flow system..." eqs = [ - connect(mass1.hp, th_resistor.hp1, th_conductor.hp1) - connect(th_conductor.hp2, flow_src.hp, hf_sensor1.hp1, hf_sensor2.hp1) - connect(th_resistor.hp2, hf_sensor1.hp2, hf_sensor2.hp2, th_ground.hp) + connect(mass1.port, th_resistor.port_a, th_conductor.port_a) + connect(th_conductor.port_b, flow_src.port, hf_sensor1.port_a, hf_sensor2.port_a) + connect(th_resistor.port_b, hf_sensor1.port_b, hf_sensor2.port_b, th_ground.port) ] @named h2 = ODESystem(eqs, t, systems=[mass1, hf_sensor1, hf_sensor2, @@ -87,13 +85,13 @@ end th_resistor.Q_flow => 1.0 ] prob = ODEProblem(sys, u0, (0, 3.0)) - sol = solve(prob, Rosenbrock23()) + sol = solve(prob, Rodas4()) - @test sol[th_conductor.T]*G == sol[th_conductor.Q_flow] - @test sol[th_conductor.Q_flow] ≈ sol[hf_sensor1.Q_flow] + sol[flow_src.hp.Q_flow] + @test sol[th_conductor.dT] .* G == sol[th_conductor.Q_flow] + @test sol[th_conductor.Q_flow] ≈ sol[hf_sensor1.Q_flow] + sol[flow_src.port.Q_flow] - @test sol[mass1.T] == sol[th_resistor.T] - @test sol[th_resistor.T]./R ≈ sol[th_resistor.Q_flow] + @test sol[mass1.T] == sol[th_resistor.port_a.T] + @test sol[th_resistor.dT] ./ R ≈ sol[th_resistor.Q_flow] end @@ -110,10 +108,10 @@ end @info "Building a piston-cylinder..." eqs = [ - connect(gas_tem.hp, gas.solidport) - connect(gas.fluidport, wall.hp1) - connect(wall.hp2, coolant.fluidport) - connect(coolant.solidport, coolant_tem.hp) + connect(gas_tem.port, gas.solidport) + connect(gas.fluidport, wall.port_a) + connect(wall.port_b, coolant.fluidport) + connect(coolant.solidport, coolant_tem.port) ] @named piston = ODESystem(eqs, t, systems=[gas_tem, wall, gas, coolant, coolant_tem]) sys = structural_simplify(piston) @@ -123,7 +121,7 @@ end wall.Q_flow => 10.0 ] prob = ODEProblem(sys, u0, (0, 3.0)) - sol = solve(prob, Rosenbrock23()) + sol = solve(prob, Rodas4()) # Heat-flow-rate is equal in magnitude # and opposite in direction @@ -141,56 +139,60 @@ end @named gas_tem = FixedTemperature(T=Tᵧ) @named coolant_tem = FixedTemperature(T=Tᵪ) @named radiator = BodyRadiation(G=G) - @named ground = ThermalGround() @named dissipator = ConvectiveConductor(G=10) + @named mass = HeatCapacitor(C=10) @info "Building a radiator..." eqs = [ - connect(gas_tem.hp, radiator.hp1, base.hp1, dissipator.solidport) - connect(base.hp2, radiator.hp2, coolant_tem.hp, dissipator.fluidport) + connect(gas_tem.port, radiator.port_a, base.port_a, dissipator.solid, mass.port) + connect(coolant_tem.port, base.port_b, radiator.port_b, dissipator.fluid) ] - @named rad = ODESystem(eqs, t, systems=[base, gas_tem, radiator, dissipator, coolant_tem]) + @named rad = ODESystem(eqs, t, systems=[base, gas_tem, radiator, dissipator, coolant_tem, mass]) sys = structural_simplify(rad) u0 = [ base.Q_flow => 10 dissipator.Q_flow => 10 + mass.T => Tᵧ ] prob = ODEProblem(sys, u0, (0, 3.0)) - sol = solve(prob, Rosenbrock23()) + sol = solve(prob, Rodas4()) - @test sol[dissipator.dT] == sol[radiator.hp1.T] - sol[radiator.hp2.T] + @test sol[dissipator.dT] == sol[radiator.port_a.T] - sol[radiator.port_b.T] rad_Q_flow = G*σ*(Tᵧ^4 - Tᵪ^4) @test sol[radiator.Q_flow] == fill(rad_Q_flow, length(sol[radiator.Q_flow])) end @testset "Thermal Collector" begin - @named flow_src = FixedHeatFlow(Q_flow=50, α=100) + @named flow_src = FixedHeatFlow(Q_flow=50, alpha=100) @named hf_sensor = HeatFlowSensor() - @named th_ground = ThermalGround() - @named collector = ThermalCollector(N=2) + @named th_ground = FixedTemperature(T=0) + @named collector = ThermalCollector(m=2) @named th_resistor = ThermalResistor(R=10) @named tem_src = FixedTemperature(T=10) + @named mass = HeatCapacitor(C=10) @info "Building a heat collector..." eqs = [ - connect(flow_src.hp, collector.hp1, th_resistor.hp1) - connect(tem_src.hp, collector.hp2) - connect(hf_sensor.hp1, collector.collector_port) - connect(hf_sensor.hp2, th_ground.hp, th_resistor.hp2) + connect(flow_src.port, collector.port_a1, th_resistor.port_a) + connect(tem_src.port, collector.port_a2) + connect(hf_sensor.port_a, collector.port_b) + connect(hf_sensor.port_b, mass.port, th_resistor.port_b) + connect(mass.port, th_ground.port) ] @named coll = ODESystem(eqs, t, systems=[hf_sensor,flow_src, tem_src, - collector, th_resistor]) + collector, th_resistor, mass]) sys = structural_simplify(coll) u0 = [ - th_resistor.Q_flow => 1.0 + th_resistor.Q_flow => 1.0, + mass.T => 0.0, ] prob = ODEProblem(sys, u0, (0, 3.0)) - sol = solve(prob, Rosenbrock23()) + sol = solve(prob, Rodas4()) - @test sol[collector.collector_port.Q_flow] + sol[collector.hp1.Q_flow] + sol[collector.hp2.Q_flow] == - zeros(length(sol[collector.collector_port.Q_flow])) - @test sol[collector.collector_port.T] == sol[collector.hp1.T] == sol[collector.hp2.T] + @test sol[collector.port_b.Q_flow] + sol[collector.port_a1.Q_flow] + sol[collector.port_a2.Q_flow] == + zeros(length(sol[collector.port_b.Q_flow])) + @test sol[collector.port_b.T] == sol[collector.port_a1.T] == sol[collector.port_a2.T] end \ No newline at end of file diff --git a/test/analog.jl b/test/analog.jl deleted file mode 100644 index b6f2a0b7a..000000000 --- a/test/analog.jl +++ /dev/null @@ -1,293 +0,0 @@ -using ModelingToolkitStandardLibrary.Electrical, ModelingToolkit, OrdinaryDiffEq, Test -using ModelingToolkitStandardLibrary.Electrical: _step, _square_wave, _triangular_wave, - _cos_wave, _damped_sine_wave, _ramp -# using Plots - -@parameters t -@named ground = Ground() -R = 5 -@named resistor = Resistor(R=R) - -@info "Testing the voltage sources..." -@testset "voltage sources" begin - @named resistor = Resistor(R=R) - @named voltage_sensor = VoltageSensor() - offset = 1 - freq = 50 - @named sinesource = SineVoltage(offset=offset, amplitude=1.0, frequency=freq, starttime=.5, phase=0.0) - @named constsource = ConstantVoltage(V=1.0) - @named stepsource = StepVoltage(height=10, offset=1.0, starttime=1.0) - sources = [constsource, sinesource, stepsource] - for source in sources - rc_eqs = [ - connect(voltage_sensor.p, source.p, resistor.p) - connect(voltage_sensor.n, source.n, resistor.n, ground.g) - ] - @named rc_model = ODESystem(rc_eqs, t, systems = [resistor, source, voltage_sensor, ground]) - sys = structural_simplify(rc_model) - u0 = [ - resistor.p.i => 10.0 - constsource.v => 1.0 - ] - prob = ODEProblem(sys, u0, (0, 2.0)) - sol = solve(prob, Rosenbrock23()) - - if source == constsource - @test sol[voltage_sensor.v][1:20000:end] ≈ ones(length(sol.t[1:20000:end])) atol=1e-3 - elseif source == stepsource - @test sol[voltage_sensor.v][1:20000:end] ≈ [(t < 1) ? offset : offset + height for t in sol.t[1:20000:end]] atol=1e-3 - else - @test sol[voltage_sensor.v][1:20000:end] ≈ [(t < 0.5) ? offset : offset + sin(2π*freq*(t .- 0.5)) for t in sol.t[1:20000:end]] atol=1e-6 - end - end -end - -@info "Testing the sensors..." -@testset "sensors" begin - offset = 1 - freq = 50 - @named source = SineVoltage(offset=offset, amplitude=1.0, frequency=freq, starttime=.5, phase=0.0) - - @testset "current" begin - @named current_sensor = CurrentSensor() - - rc_eqs = [ - connect(source.p, resistor.p) - connect(resistor.n, current_sensor.p) - connect(current_sensor.n, source.n, ground.g) - ] - @named rc_model = ODESystem(rc_eqs, t, systems = [resistor, source, current_sensor, ground]) - sys = structural_simplify(rc_model) - u0 = [ - resistor.p.i => 10.0 - ] - prob = ODEProblem(sys, u0, (0, 2.0)) - sol = solve(prob, Rosenbrock23()) - - @test sol[current_sensor.i][1:20000:end] ≈ [(t < 0.5) ? offset/R : (offset + sin(2π*freq*(t .- 0.5)))/R for t in sol.t[1:20000:end]] atol=1e-6 - end - - @testset "potential" begin - @named potential_sensor1 = PotentialSensor() - @named potential_sensor2 = PotentialSensor() - - rc_eqs = [ - connect(source.p, resistor.p, potential_sensor1.p) - connect(resistor.n, potential_sensor2.p, source.n, ground.g) - ] - @named rc_model = ODESystem(rc_eqs, t, systems = [resistor, source, potential_sensor1, potential_sensor2, ground]) - sys = structural_simplify(rc_model) - u0 = [ - resistor.p.i => 10.0 - ] - prob = ODEProblem(sys, u0, (0, 2.0)) - sol = solve(prob, Rosenbrock23()) - @test sol[potential_sensor1.phi][1:20000:end] ≈ [(t < 0.5) ? offset : offset + sin(2π*freq*(t .- 0.5)) for t in sol.t[1:20000:end]] atol=1e-6 - @test iszero(sol[potential_sensor2.phi][1:20000:end]) - end - - @testset "voltage" begin - @named voltage_sensor = VoltageSensor() - - rc_eqs = [ - connect(source.p, resistor.p, voltage_sensor.p) - connect(voltage_sensor.n, source.n, resistor.n, ground.g) - ] - @named rc_model = ODESystem(rc_eqs, t, systems = [resistor, source, voltage_sensor, ground]) - sys = structural_simplify(rc_model) - u0 = [ - resistor.p.i => 10.0 - ] - prob = ODEProblem(sys, u0, (0, 2.0)) - sol = solve(prob, Rosenbrock23()) - @test sol[voltage_sensor.v][1:20000:end] ≈ [(t < 0.5) ? offset : offset + sin(2π*freq*(t .- 0.5)) for t in sol.t[1:20000:end]] atol=1e-6 - end - - @testset "power" begin - @named power_sensor = PowerSensor() - - rc_eqs = [ - connect(source.p, resistor.p, power_sensor.pv) - connect(power_sensor.nv, resistor.n, power_sensor.nc) - connect(power_sensor.pc, source.n, ground.g) - ] - @named rc_model = ODESystem(rc_eqs, t, systems = [resistor, source, power_sensor, ground]) - sys = structural_simplify(rc_model) - u0 = [ - resistor.p.i => 10.0 - ] - prob = ODEProblem(sys, u0, (0, 2.0)) - sol = solve(prob, Rosenbrock23()) - @test sol[power_sensor.power][1:20000:end] ≈ [(t < 0.5) ? offset^2/R : (offset + sin(2π*freq*(t .- 0.5)))^2 / R for t in sol.t[1:20000:end]] atol=1e-6 - end - - @testset "multi" begin - @named multi_sensor = MultiSensor() - - rc_eqs = [ - connect(source.p, resistor.p, multi_sensor.pv) - connect(multi_sensor.nv, resistor.n, multi_sensor.nc) - connect(multi_sensor.pc, source.n, ground.g) - ] - @named rc_model = ODESystem(rc_eqs, t, systems = [resistor, source, multi_sensor, ground]) - sys = structural_simplify(rc_model) - u0 = [ - resistor.p.i => 10.0 - ] - prob = ODEProblem(sys, u0, (0, 2.0)) - sol = solve(prob, Rosenbrock23()) - @test sol[multi_sensor.i][1:20000:end] ≈ [(t < 0.5) ? offset/R : (offset + sin(2π*freq*(t .- 0.5)))/R for t in sol.t[1:20000:end]] atol=1e-6 - @test sol[multi_sensor.v][1:20000:end] ≈ [(t < 0.5) ? offset : offset + sin(2π*freq*(t .- 0.5)) for t in sol.t[1:20000:end]] atol=1e-6 - end -end - -@info "Testing the inductor..." -@testset "Inductor" begin - freq, offset = 5000, 1 - @named sinesource = SineVoltage(offset=offset, starttime=0.5, amplitude=1.0, frequency=5, phase=0.0) - @named l1 = Inductor() - @named l2 = Inductor() - @named voltage_sensor = VoltageSensor() - l_eqs = [ - connect(voltage_sensor.p, sinesource.p, l1.p) - connect(l1.n, l2.p) - connect(voltage_sensor.n, sinesource.n, l2.n, ground.g) - ] - @named l_model = ODESystem(l_eqs, t, systems = [l1, l2, sinesource, voltage_sensor, ground]) - sys = structural_simplify(l_model) - equations(sys) - u0 = [ - l1.p.i => 10.0 - sinesource.v => 1.0 - l2.v => 0.0 - ] - prob = ODEProblem(sys, u0, (0, 10.0)) - sol = solve(prob, Rosenbrock23()) - - @test sol[l1.v] + sol[l2.v] ≈ sol[voltage_sensor.v] -end - -@info "Constructing an integrator circuit..." -# This tests Capacitor, IdealOpAmp -@testset "Integrator" begin - @named ground = Ground() - @named res1 = Resistor() - @named c1 = Capacitor() - @named opamp = IdealOpAmp() - @named square = SquareVoltage() - - in_eqs = [ - connect(square.p, res1.p) - connect(res1.n, c1.p, opamp.p1) - connect(opamp.n2, c1.n) - connect(opamp.n1, ground.g, opamp.p2, square.n) - ] - @named in_model = ODESystem(in_eqs, t, systems = [res1, opamp, square, c1, ground]) - sys = structural_simplify(in_model) - u0 = [ - c1.v => 1 - res1.v => 1 - ] - prob = ODEProblem(sys, u0, (0, 10.0)) - sol = solve(prob, Rosenbrock23()) - @test sol[opamp.v2] == sol[c1.v] # Not a great one however. Rely on the plot - - # plt = plot(sol) - # savefig(plt, "integrator") -end - -@info "Testing voltage function generators..." -@testset "Voltage function generators" begin - st, o, h, f, A, et, ϕ, d, δ = 0.7, 1.25, 3, 2, 2.5, 2.5, π/4, 0.1, 0.0001 - - @named ground = Ground() - @named res = Resistor() - @named voltage_sensor = VoltageSensor() - @named vstep = StepVoltage(starttime=st, offset=o, height=h) - @named vsquare = SquareVoltage(offset=o, starttime=st, amplitude=A, frequency=f) - @named vtri = TriangularVoltage(offset=o, starttime=st, amplitude=A, frequency=f) - # @named vsawtooth = SawToothVoltage(amplitude=A, starttime=st, frequency=f, offset=o) - @named vcosine = CosineVoltage(offset=o, amplitude=A, frequency=f, starttime=st, phase=ϕ) - @named vdamped_sine = DampedSineVoltage(offset=o, amplitude=A, frequency=f, starttime=st, phase=ϕ, damping_coef=d) - @named vramp = RampVoltage(offset=o, starttime=st, endtime=et, height=h) - - vsources = [vtri, vsquare, vstep, vcosine, vdamped_sine, vramp] - waveforms(i, x) = getindex([o .+ (x .> st) .* _triangular_wave.(x, δ, f, A, st), - o .+ (x .> st) .* _square_wave.(x, δ, f, A, st), - o .+ _step.(x, δ, h, st), - # o .+ (x .> st). * _sawtooth_wave.(x, δ, f, A, st), - o .+ (x .> st) .* _cos_wave.(x, f, A, st, ϕ), - o .+ (x .> st) .* _damped_sine_wave.(x, f, A, st, ϕ, d), - o .+ _ramp.(x, δ, st, et, h)], i) - for i in 1:length(vsources) - vsource = vsources[i] - @info Symbolics.getname(vsource) - eqs = [ - connect(vsource.p, voltage_sensor.p, res.p) - connect(ground.g, voltage_sensor.n, vsource.n, res.n) - ] - @named vmodel = ODESystem(eqs, t, systems = [voltage_sensor, res, vsource, ground]) - vsys = structural_simplify(vmodel) - - u0 = [ - vsource.v => 1 - res.v => 1 - ] - - prob = ODEProblem(vsys, u0, (0, 10.0)) - vsol = solve(prob, dt=0.1, Rosenbrock23()) - - @test vsol[vsource.v][1150:end] ≈ waveforms(i, vsol.t)[1150:end] atol=1e-1 - # For visual inspection - # plt = plot(sol) - # savefig(plt, "test_current_$(Symbolics.getname(source))") - end -end - -@info "Testing the Current generators..." -@testset "Current function generators" begin - st, o, h, f, A, et, ϕ, d, δ = 0.7, 1.25, 3, 2, 2.5, 2.5, π/4, 0.1, 0.0001 - - @named ground = Ground() - @named res = Resistor() - @named current_sensor = CurrentSensor() - @named istep = StepCurrent(starttime=st, offset=o, height=h) - @named isquare = SquareCurrent(offset=o, starttime=st, amplitude=A, frequency=f) - @named itri = TriangularCurrent(offset=o, starttime=st, amplitude=A, frequency=f) - # @named isawtooth = SawToothCurrent(amplitude=A, starttime=st, frequency=f, offset=o) - @named icosine = CosineCurrent(offset=o, amplitude=A, frequency=f, starttime=st, phase=ϕ) - @named idamped_sine = DampedSineCurrent(offset=o, amplitude=A, frequency=f, starttime=st, phase=ϕ, damping_coef=d) - @named iramp = RampCurrent(offset=o, starttime=st, endtime=et, height=h) - - isources = [itri, isquare, istep, icosine, idamped_sine, iramp] - waveforms(i, x) = getindex([o .+ (x .> st) .* _triangular_wave.(x, δ, f, A, st), - o .+ (x .> st) .* _square_wave.(x, δ, f, A, st), - o .+ _step.(x, δ, h, st), - # o .+ (x .> st). * _sawtooth_wave.(x, δ, f, A, st), - o .+ (x .> st) .* _cos_wave.(x, f, A, st, ϕ), - o .+ (x .> st) .* _damped_sine_wave.(x, f, A, st, ϕ, d), - o .+ _ramp.(x, δ, st, et, h)], i) - - for i in 1:length(isources) - isource = isources[i] - eqs = [ - connect(isource.p, current_sensor.n) - connect(current_sensor.p, res.n) - connect(isource.n, res.p) - ] - @named model = ODESystem(eqs, t, systems = [current_sensor, isource, res]) - isys = alias_elimination(model) - - u0 = [ - isource.i => 1 - res.v => 1 - ] - prob = ODEProblem(isys, u0, (0, 2.0)) - sol = solve(prob, Rosenbrock23()) - - @test sol[isource.i][1150:end] ≈ waveforms(i, sol.t)[1150:end] atol=1e-1 - # For visual inspection - # plt = plot(sol) - # savefig(plt, "test_current_$(Symbolics.getname(source))") - end -end \ No newline at end of file diff --git a/test/demo.jl b/test/demo.jl index f9cca3a37..42b7c0856 100644 --- a/test/demo.jl +++ b/test/demo.jl @@ -1,25 +1,27 @@ using ModelingToolkitStandardLibrary.Electrical, ModelingToolkit, OrdinaryDiffEq #, Plots +using Test -R = 1.0 -C = 1.0 -V = 1.0 -@named resistor = Resistor(R=R) -@named capacitor = Capacitor(C=C) -@named source = ConstantVoltage(V=V) -@named ground = Ground() +@testset "RC demo" begin + R = 1.0 + C = 1.0 + V = 1.0 + @parameters t + @named resistor = Resistor(R=R) + @named capacitor = Capacitor(C=C) + @named source = ConstantVoltage(V=V) + @named ground = Ground() -rc_eqs = [ - connect(source.p, resistor.p) - connect(resistor.n, capacitor.p) - connect(capacitor.n, source.n, ground.g) - ] + rc_eqs = [ + connect(source.p, resistor.p) + connect(resistor.n, capacitor.p) + connect(capacitor.n, source.n, ground.g) + ] -@named rc_model = ODESystem(rc_eqs, systems=[resistor, capacitor, source, ground]) -sys = structural_simplify(rc_model) -u0 = [ - capacitor.v => 0.0 - capacitor.p.i => 0.0 - ] -prob = ODAEProblem(sys, u0, (0, 10.0)) -sol = solve(prob, Tsit5()) -#plot(sol) + @named rc_model = ODESystem(rc_eqs, t, systems=[resistor, capacitor, source, ground]) + sys = structural_simplify(rc_model) + prob = ODAEProblem(sys, Pair[], (0, 10.0)) + sol = solve(prob, Tsit5()) + #plot(sol) + + @test isapprox(sol[capacitor.v][end], V, atol=1e-2) +end \ No newline at end of file diff --git a/test/magnetic.jl b/test/magnetic.jl deleted file mode 100644 index 7905a6557..000000000 --- a/test/magnetic.jl +++ /dev/null @@ -1,7 +0,0 @@ -using ModelingToolkitStandardLibrary.Magnetic, ModelingToolkit, OrdinaryDiffEq, Test - -@parameters t -@named ground = Ground() - -@info "Testing basic magnetic components..." - diff --git a/test/runtests.jl b/test/runtests.jl index c55d5571d..f606e70f1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,9 +1,22 @@ using SafeTestsets -@safetestset "Blocks: math" begin include("test_math.jl") end -@safetestset "Blocks: nonlinear" begin include("test_nonlinear.jl") end -@safetestset "Blocks: continuous" begin include("test_continuous.jl") end -@safetestset "Analog Circuits" begin include("analog.jl") end -#@safetestset "Digital Circuits" begin include("digital.jl") end +# Blocks +@safetestset "Blocks: math" begin include("Blocks/math.jl") end +@safetestset "Blocks: nonlinear" begin include("Blocks/nonlinear.jl") end +@safetestset "Blocks: continuous" begin include("Blocks/continuous.jl") end +@safetestset "Blocks: sources" begin include("Blocks/sources.jl") end + +# Electrical +@safetestset "Analog Circuits" begin include("Electrical/analog.jl") end +#@safetestset "Digital Circuits" begin include("Electrical/digital.jl") end @safetestset "RC Circuit Demo" begin include("demo.jl") end -@safetestset "Thermal Circuits" begin include("thermal.jl") end + +# Thermal +@safetestset "Thermal Circuits" begin include("Thermal/thermal.jl") end +@safetestset "Thermal Demo" begin include("Thermal/demo.jl") end + +# Magnetic +@safetestset "Magnetic" begin include("Magnetic/magnetic.jl") end + +# Mechanical +@safetestset "Mechanical" begin include("Mechanical/rotational.jl") end diff --git a/test/test_continuous.jl b/test/test_continuous.jl deleted file mode 100644 index 8e00b6d69..000000000 --- a/test/test_continuous.jl +++ /dev/null @@ -1,269 +0,0 @@ -using ModelingToolkit, ModelingToolkitStandardLibrary, OrdinaryDiffEq -using ModelingToolkitStandardLibrary.Blocks -using ModelingToolkitStandardLibrary.Blocks: Dₜ, t - -#= -Testing strategy: -The general strategy is to test systems using simple intputs where the solution -is known on closed form. For algebraic systems (without differential variables), -an integrator with a constant input is often used together with the system under test. -=# - -@testset "Constant" begin - @info "Testing Constant" - @named c = Constant(42) - sys = structural_simplify(c) - prob = ODEProblem(sys, Pair[], (0.0, 1.0)) - @test_broken begin # StackOverflowError: - sol = solve(prob, Rosenbrock23()) - @test all(==(42), sol[c.y]) - end - @named int = Integrator(; k=1) - @named iosys = ODESystem([int.u~c.y], t, systems=[c, int]) - sys = structural_simplify(iosys) - prob = ODEProblem(sys, Pair[], (0.0, 1.0)) - sol = solve(prob, Rosenbrock23(), saveat=0:0.1:1) - @test sol[int.y] ≈ 42 .* (0:0.1:1) -end - -@testset "Integrator" begin - @info "Testing Integrator" - for k = [0.1, 1, 10] - @named int = Integrator(; k) - @test count(ModelingToolkit.isinput, states(int)) == 1 - @test count(ModelingToolkit.isoutput, states(int)) == 1 - @named iosys = ODESystem([int.u~1], t, systems=[int]) - sys = structural_simplify(iosys) - stateind(sym) = findfirst(isequal(sym),states(sys)) - prob = ODEProblem(sys, Pair[int.u=>1.], (0.0, 1.0)) - sol = solve(prob, Rosenbrock23(), saveat=0:0.1:1) - @test sol[int.y] ≈ k .* (0:0.1:1) - end -end - -@testset "Derivative" begin - @info "Testing Derivative" - - #= Derivative - The test output below is generated by - using ControlSystems - sys = ss(-1/T, 1/T, -k/T, k/T) - tv = 0:0.5:10 - u = (x,t)->[sin(t)] - y = vec(lsim(sys, u, tv, alg=Rosenbrock23())[1]) - =# - k = 1 - T = 0.1 - - y01 = [0.0, 0.9096604391560481, 0.6179369162956885, 0.16723968919320775, -0.3239425882305049, -0.7344654437585882, -0.9662915429884467, -0.9619031643363591, -0.7219189996926385, -0.3046471954239838, 0.18896274787342904, 0.6325612488150467, 0.923147635361496, 0.9882186461533009, 0.8113758856575801, 0.4355269842556595, -0.05054266121798534, -0.5180957852231662, -0.8615644854197235, -0.994752654263345, -0.8845724777509947] - - y1 = [0.0, 0.37523930001382705, 0.5069379343173124, 0.422447016206449, 0.17842742193310424, -0.14287580928455357, -0.44972307981519677, -0.6589741190943343, -0.7145299845867902, -0.5997749247850142, -0.34070236779586216, -5.95731929625698e-5, 0.33950710748637825, 0.595360048429, 0.7051403889991136, 0.6421181090255983, 0.4214753349401378, 0.09771852881756515, -0.24995564964733635, -0.5364893060362096, -0.6917461951831227] - - y10 = [0.0, 0.04673868865158038, 0.07970450452536708, 0.09093906605247397, 0.07779607227750623, 0.04360203242101193, -0.0031749143460660587, -0.050989771426848074, -0.08804727520541561, -0.10519046453331109, -0.09814083278784949, -0.06855209962041757, -0.023592611490189652, 0.025798926487949535, 0.0675952553752348, 0.0916256775597053, 0.09206230764744555, 0.06885879535935949, 0.027748930190142837, -0.021151336671582116, -0.06582115823326284] - - for k = [0.1, 1, 10], (T,y) = zip([0.1, 1, 10], [y01, y1, y10]) - @named der = Derivative(; k, T) - @test count(ModelingToolkit.isinput, states(der)) == 1 - @test count(ModelingToolkit.isoutput, states(der)) == 1 - @named iosys = ODESystem([der.u~sin(t)], t, systems=[der]) - sys = structural_simplify(iosys) - stateind(sym) = findfirst(isequal(sym),states(sys)) - prob = ODEProblem(sys, Pair[der.u=>0., der.x=>0], (0.0, 10.0)) - sol = solve(prob, Rosenbrock23(), saveat=0:0.5:10) - # plot([sol[der.y] k.*y]) |> display - @test sol[der.y] ≈ k .* y rtol=1e-2 - end -end - -@testset "FirstOrder" begin - @info "Testing FirstOrder" - - for k = [0.1, 1, 10], T = [0.1, 1, 10] - @named fo = FirstOrder(; k, T) - @test count(ModelingToolkit.isinput, states(fo)) == 1 - @test count(ModelingToolkit.isoutput, states(fo)) == 1 - @named iosys = ODESystem([fo.u~1], t, systems=[fo]) - sys = structural_simplify(iosys) - prob = ODEProblem(sys, Pair[fo.u=>1., fo.x=>0], (0.0, 10.0)) - sol = solve(prob, Rosenbrock23(), saveat=0:0.1:10) - y = k .* (1 .- exp.(.-sol.t ./ T)) # Known solution to first-order system - # plot([sol[fo.y] y]) |> display - @test sol[fo.y] ≈ y rtol=1e-3 - end - -end - -@testset "SecondOrder" begin - @info "Testing SecondOrder" - - # The impulse response of a second-order system with damping d follows the equations below - function so(t,w,d) - val = if d == 0 - 1/w * sin(w*t) - elseif d < 1 - 1/(w*sqrt(1-d^2)) * exp(-d*w*t) * sin(w*sqrt(1-d^2)*t) - elseif d == 1 - t*exp(-w*t) - else - 1/(w*sqrt(d^2-1)) * exp(-d*w*t) * sinh(w*sqrt(d^2-1)*t) - end - val - end - - w = 1 - d = 0.5 - for k = [0.1, 1, 10], w = [0.1, 1, 10], d = [0, 0.01, 0.1, 1, 1.1] - @named sos = SecondOrder(; k, w, d) - @test count(ModelingToolkit.isinput, states(sos)) == 1 - @test count(ModelingToolkit.isoutput, states(sos)) == 1 - @named iosys = ODESystem([sos.u~0], t, systems=[sos]) - sys = structural_simplify(iosys) - prob = ODEProblem(sys, Pair[sos.u=>0.0, sos.xd=>1.0], (0.0, 10.0)) # set initial derivative state to 1 to simulate an impulse response - sol = solve(prob, Rosenbrock23(), saveat=0:0.1:10, reltol=1e-6) - y = so.(sol.t,w,d)# Known solution to second-order system - # plot([sol[sos.y] y]) |> display - @test sum(abs2, sol[sos.y] - y) < 1e-4 - end -end - -@testset "PID" begin - @info "Testing PID" - - k = 2 - Ti = 0.5 - Td = 0.7 - wp = 1 - wd = 1 - Ni = √(Td / Ti) - Nd = 12 - y_max = Inf - y_min = -Inf - u_r = sin(t) - u_y = 0 - function solve_with_input(; u_r, u_y, - controller = PID(; k, Ti, Td, wp, wd, Ni, Nd, y_max, y_min, name=:controller) - ) - @test count(ModelingToolkit.isinput, states(controller)) == 5 # 2 in PID, 1 sat, 1 I, 1 D - @test count(ModelingToolkit.isoutput, states(controller)) == 4 - # TODO: check number of unbound inputs when available, should be 2 - @named iosys = ODESystem([controller.u_r~u_r, controller.u_y~u_y], t, systems=[controller]) - sys = structural_simplify(iosys) - prob = ODEProblem(sys, Pair[], (0.0, 10.0)) - sol = solve(prob, Rosenbrock23(), saveat=0:0.1:10) - controller, sys, sol - end - - # linearity in u_r - controller, sys, sol1 = solve_with_input(u_r=sin(t), u_y=0) - controller, sys, sol2 = solve_with_input(u_r=2sin(t), u_y=0) - @test sum(abs, sol1[controller.ea]) < eps() # This is the acutator model error due to saturation - @test 2sol1[controller.y] ≈ sol2[controller.y] rtol=1e-3 # linearity in u_r - - # linearity in u_y - controller, sys, sol1 = solve_with_input(u_y=sin(t), u_r=0) - controller, sys, sol2 = solve_with_input(u_y=2sin(t), u_r=0) - @test sum(abs, sol1[controller.ea]) < eps() # This is the acutator model error due to saturation - @test 2sol1[controller.y] ≈ sol2[controller.y] rtol=1e-3 # linearity in u_y - - # zero error - controller, sys, sol1 = solve_with_input(u_y=sin(t), u_r=sin(t)) - @test sum(abs, sol1[controller.y]) ≈ 0 atol=sqrt(eps()) - - # test saturation - controller, sys, sol1 = solve_with_input(; u_r=10sin(t), u_y=0, - controller = PID(; k, Ti, Td, wp, wd=0, Ni, Nd, y_max=10, y_min=-10, name=:controller) - ) - @test extrema(sol1[controller.y]) == (-10, 10) - - - # test P set-point weighting - controller, sys, sol1 = solve_with_input(; u_r=sin(t), u_y=0, - controller = PID(; k, Ti, Td, wp=0, wd, Ni, Nd, y_max, y_min, name=:controller) - ) - @test sum(abs, sol1[controller.ep]) ≈ 0 atol=sqrt(eps()) - - # test D set-point weighting - controller, sys, sol1 = solve_with_input(; u_r=sin(t), u_y=0, - controller = PID(; k, Ti, Td, wp, wd=0, Ni, Nd, y_max, y_min, name=:controller) - ) - @test sum(abs, sol1[controller.ed]) ≈ 0 atol=sqrt(eps()) - - - # zero integral gain - controller, sys, sol1 = solve_with_input(; u_r=sin(t), u_y=0, - controller = PID(; k, Ti=false, Td, wp, wd, Ni, Nd, y_max, y_min, name=:controller) - ) - @test isapprox(sum(abs, sol1[controller.I.y]), 0, atol=sqrt(eps())) - - - # zero derivative gain - @test_skip begin # During the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a non-finite number: [5] - controller, sys, sol1 = solve_with_input(; u_r=sin(t), u_y=0, - controller = PID(; k, Ti, Td=false, wp, wd, Ni, Nd, y_max, y_min, name=:controller) - ) - @test isapprox(sum(abs, sol1[controller.D.y]), 0, atol=sqrt(eps())) - end - - # Tests below can be activated when the concept of unbound_inputs exists in MTK - # @test isequal(Set(unbound_inputs(controller)), @nonamespace(Set([controller.u_r, controller.u_y]))) - # @test isempty(unbound_inputs(sys)) - # @test isequal(bound_inputs(sys), inputs(sys)) - # @test isequal( - # Set(bound_inputs(sys)), - # Set([controller.u_r, controller.u_y, controller.I.u, controller.D.u, controller.sat.u]) - # ) -end - -## Additional test of PID controller using ControlSystems -# using ControlSystems -# kd = 1 -# Nd = 12 -# Td = 1 -# T = Td/Nd -# Cd = ss(-1/T, 1/T, -kd/T, kd/T) |> tf - -# C = ControlSystems.pid(kp=10, ki=1, kd=0, series=true, time=true) + 10*Cd -# P = tf(1,[1, 0])^2 -# L = ss(P*C) - -# @named controller = PID(k=10, Ti=1, Td=1) -# @named plant = Blocks.StateSpace(ssdata(ss(P))...) -# @named iosys = ODESystem([ -# controller.u_r~1, -# controller.u_y~plant.y[1], -# controller.y~plant.u[1] -# ], t, systems=[controller, plant]) -# sys = structural_simplify(iosys) -# prob = ODEProblem(sys, Pair[], (0.0, 6)) -# sol = solve(prob, Rosenbrock23()) - -# res = step(feedback(L), sol.t) -# y = res.y[:] -# plot(res) -# plot!(sol, vars=[plant.y[1]]) -# @test sol[plant.y[1]] ≈ y rtol = 1e-3 -## - -@testset "StateSpace" begin - @info "Testing StateSpace" - - A = [0 1; 0 0] - B = [0, 1] - C = [1 0] - D = 0 - @named sys = Blocks.StateSpace(A,B,C,D) - @test count(ModelingToolkit.isinput, states(sys)) == 1 - @test count(ModelingToolkit.isoutput, states(sys)) == 1 - @named iosys = ODESystem([sys.u[1] ~ 1], t, systems=[sys]) - iosys = structural_simplify(iosys) - prob = ODEProblem(iosys, Pair[], (0.0, 1.0)) - sol = solve(prob, Rosenbrock23(), saveat=0:0.1:1) - @test sol[sys.x[2]] ≈ (0:0.1:1) - @test sol[sys.x[1]] ≈ sol[sys.y[1]] - - - D = randn(2, 2) # If there's only a `D` matrix, the result is a matrix gain - @named sys = Blocks.StateSpace([],[],[],D) - gain = Blocks.Gain(D, name=:sys) - @test sys == gain -end \ No newline at end of file diff --git a/test/test_math.jl b/test/test_math.jl deleted file mode 100644 index f54a550e0..000000000 --- a/test/test_math.jl +++ /dev/null @@ -1,64 +0,0 @@ -using ModelingToolkit, OrdinaryDiffEq -using ModelingToolkitStandardLibrary.Blocks -using ModelingToolkitStandardLibrary.Blocks: Dₜ, t - -#= -Testing strategy: -The general strategy is to test systems using simple intputs where the solution is known on closed form. For algebraic systems (without differential variables), an integrator with a constant input is often used together with the system under test. -=# - -@testset "Gain" begin - @info "Testing Gain" - @named c = Gain(42) - @named int = Integrator(; k=1) - @named iosys = ODESystem([int.u~c.y, c.u~1], t, systems=[c, int]) - sys = structural_simplify(iosys) - prob = ODEProblem(sys, Pair[], (0.0, 1.0)) - sol = solve(prob, Rosenbrock23(), saveat=0:0.1:1) - @test sol[int.y] ≈ 42 .* (0:0.1:1) - - # Matrix gain - @named c = Gain([2 0; 0 3]) - ints = [Integrator(; k=1, name=Symbol("int$i")) for i in 1:2] - @named iosys = ODESystem([ - ints[1].u~c.y[1], - ints[2].u~c.y[2], - c.u[1]~1, - c.u[2]~2, - ], t, systems=[c; ints]) - sys = structural_simplify(iosys) - prob = ODEProblem(sys, Pair[], (0.0, 1.0)) - sol = solve(prob, Rosenbrock23(), saveat=0:0.1:1) - @test sol[ints[1].y] ≈ 2 .* (0:0.1:1) # 2 * 1 - @test sol[ints[2].y] ≈ 6 .* (0:0.1:1) # 3 * 2 -end - - -@testset "Sum" begin - @info "Testing Sum" - @named s = Sum(2) - ints = [Integrator(; k=1, name=Symbol("int$i")) for i in 1:2] - @named iosys = ODESystem([ - ints[1].u~1, - ints[2].u~2, - ints[1].y~s.u[1], - ints[2].y~s.u[2], - ], t, systems=[s; ints]) - sys = structural_simplify(iosys) - prob = ODEProblem(sys, Pair[], (0.0, 1.0)) - sol = solve(prob, Rosenbrock23(), saveat=0:0.1:1) - @test sol[s.y] ≈ 3 .* (0:0.1:1) - - @named s = Sum([1, -2]) - ints = [Integrator(; k=1, name=Symbol("int$i")) for i in 1:2] - @named iosys = ODESystem([ - ints[1].u~1, - ints[2].u~1, - ints[1].y~s.u[1], - ints[2].y~s.u[2], - ], t, systems=[s; ints]) - sys = structural_simplify(iosys) - prob = ODEProblem(sys, Pair[], (0.0, 1.0)) - sol = solve(prob, Rosenbrock23(), saveat=0:0.1:1) - @test sol[s.y] ≈ (1 + (-2)) .* (0:0.1:1) -end \ No newline at end of file diff --git a/test/test_nonlinear.jl b/test/test_nonlinear.jl deleted file mode 100644 index f8c41cf60..000000000 --- a/test/test_nonlinear.jl +++ /dev/null @@ -1,57 +0,0 @@ -using ModelingToolkit, ModelingToolkitStandardLibrary, OrdinaryDiffEq -using ModelingToolkitStandardLibrary.Blocks: t, Saturation, DeadZone, Integrator - -#= -Testing strategy: -The general strategy is to test systems using simple intputs where the solution is known on closed form. For algebraic systems (without differential variables), an integrator with a constant input is often used together with the system under test. -=# - -## Saturation -@testset "Saturation" begin - @info "Testing Saturation" - y_max = 0.8 - y_min = -0.6 - @named sat = Saturation(; y_min, y_max) - @test count(ModelingToolkit.isinput, states(sat)) == 1 - @test count(ModelingToolkit.isoutput, states(sat)) == 1 - @named iosys = ODESystem([sat.u~sin(t)], t, systems=[sat]) - - sys = structural_simplify(iosys) - stateind(sym) = findfirst(isequal(sym),states(sys)) - prob = ODEProblem(sys, Pair[sat.u=>0], (0.0, 2pi)) - sol = solve(prob, Rosenbrock23(), saveat=0:0.1:2pi) - y = clamp.(sin.(sol.t), y_min, y_max) - # plot([sol[sat.y] y]) - @test sol[sat.y] ≈ y rtol = 1e-6 -end - - -@testset "DeadZone" begin - @info "Testing DeadZone" - - ie = ifelse - deadzone(u, u_min, u_max) = ie(u > u_max, u-u_max, ie( u < u_min, u-u_min, 0)) - - u_max = 1 - u_min = -2 - @named dz = DeadZone(; u_min, u_max) - @named int = Integrator() - @test count(ModelingToolkit.isinput, states(dz)) == 1 - @test count(ModelingToolkit.isoutput, states(dz)) == 1 - @named iosys = ODESystem([ - int.u ~ 1 - int.y ~ dz.u - ], t, systems=[dz, int] - ) - - sys = structural_simplify(iosys) - prob = ODEProblem(sys, Pair[int.x => -3], (0.0, 5)) - sol = solve(prob, Rosenbrock23(), saveat=0:0.1:5) - y = deadzone.(sol[int.y], u_min, u_max) - @test y[1] == -3 - u_min - @test y[end] ≈ 2 - u_max - @test y[end÷2] == 0 - - # plot([y sol[dz.y] sol[dz.u]]) - @test sol[dz.y] ≈ y rtol = 1e-6 -end \ No newline at end of file