diff --git a/.github/dependabot.yml b/.github/dependabot.yml
index 700707ced..ec3b005a0 100644
--- a/.github/dependabot.yml
+++ b/.github/dependabot.yml
@@ -5,3 +5,6 @@ updates:
directory: "/" # Location of package manifests
schedule:
interval: "weekly"
+ ignore:
+ - dependency-name: "crate-ci/typos"
+ update-types: ["version-update:semver-patch", "version-update:semver-minor"]
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 92b7ab8d7..687e5f669 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -3,9 +3,13 @@ on:
pull_request:
branches:
- main
+ paths-ignore:
+ - 'docs/**'
push:
branches:
- main
+ paths-ignore:
+ - 'docs/**'
jobs:
test:
runs-on: ubuntu-latest
@@ -15,13 +19,12 @@ jobs:
- All
version:
- '1'
- - '1.6'
steps:
- - uses: actions/checkout@v3
+ - uses: actions/checkout@v4
- uses: julia-actions/setup-julia@v1
with:
version: ${{ matrix.version }}
- - uses: actions/cache@v3
+ - uses: actions/cache@v4
env:
cache-name: cache-artifacts
with:
@@ -34,6 +37,8 @@ jobs:
- uses: julia-actions/julia-buildpkg@v1
- uses: julia-actions/julia-runtest@v1
- uses: julia-actions/julia-processcoverage@v1
- - uses: codecov/codecov-action@v3
+ - uses: codecov/codecov-action@v4
with:
file: lcov.info
+ token: ${{ secrets.CODECOV_TOKEN }}
+ fail_ci_if_error: true
diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml
index bca716089..135b571b3 100644
--- a/.github/workflows/Documentation.yml
+++ b/.github/workflows/Documentation.yml
@@ -11,7 +11,7 @@ jobs:
build:
runs-on: ubuntu-latest
steps:
- - uses: actions/checkout@v3
+ - uses: actions/checkout@v4
- uses: julia-actions/setup-julia@latest
with:
version: '1'
@@ -23,6 +23,6 @@ jobs:
DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
run: julia --project=docs/ --code-coverage=user docs/make.jl
- uses: julia-actions/julia-processcoverage@v1
- - uses: codecov/codecov-action@v3
+ - uses: codecov/codecov-action@v4
with:
file: lcov.info
diff --git a/.github/workflows/Downgrade.yml b/.github/workflows/Downgrade.yml
new file mode 100644
index 000000000..93974846d
--- /dev/null
+++ b/.github/workflows/Downgrade.yml
@@ -0,0 +1,29 @@
+name: Downgrade
+on:
+ pull_request:
+ branches:
+ - main
+ paths-ignore:
+ - 'docs/**'
+ push:
+ branches:
+ - main
+ paths-ignore:
+ - 'docs/**'
+jobs:
+ test:
+ runs-on: ubuntu-latest
+ strategy:
+ matrix:
+ version: ['1']
+ steps:
+ - uses: actions/checkout@v4
+ - uses: julia-actions/setup-julia@v1
+ with:
+ version: ${{ matrix.version }}
+ - uses: julia-actions/julia-downgrade-compat@v1
+# if: ${{ matrix.version == '1.6' }}
+ with:
+ skip: Pkg,TOML
+ - uses: julia-actions/julia-buildpkg@v1
+ - uses: julia-actions/julia-runtest@v1
diff --git a/.github/workflows/Downstream.yml b/.github/workflows/Downstream.yml
index 5528e2cb0..e0ddd30b5 100644
--- a/.github/workflows/Downstream.yml
+++ b/.github/workflows/Downstream.yml
@@ -14,19 +14,19 @@ jobs:
strategy:
fail-fast: false
matrix:
- julia-version: [1,1.6]
+ julia-version: [1]
os: [ubuntu-latest]
package:
- {user: SciML, repo: ModelingToolkit.jl}
steps:
- - uses: actions/checkout@v3
+ - uses: actions/checkout@v4
- uses: julia-actions/setup-julia@v1
with:
version: ${{ matrix.julia-version }}
arch: x64
- uses: julia-actions/julia-buildpkg@latest
- name: Clone Downstream
- uses: actions/checkout@v3
+ uses: actions/checkout@v4
with:
repository: ${{ matrix.package.user }}/${{ matrix.package.repo }}
path: downstream
@@ -48,6 +48,8 @@ jobs:
exit(0) # Exit immediately, as a success
end
- uses: julia-actions/julia-processcoverage@v1
- - uses: codecov/codecov-action@v3
+ - uses: codecov/codecov-action@v4
with:
file: lcov.info
+ token: ${{ secrets.CODECOV_TOKEN }}
+ fail_ci_if_error: true
diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml
index f80d0b18b..45bd09c47 100644
--- a/.github/workflows/FormatCheck.yml
+++ b/.github/workflows/FormatCheck.yml
@@ -3,7 +3,7 @@ name: format-check
on:
push:
branches:
- - 'master'
+ - 'main'
- 'release-'
tags: '*'
pull_request:
@@ -21,7 +21,7 @@ jobs:
with:
version: ${{ matrix.julia-version }}
- - uses: actions/checkout@v3
+ - uses: actions/checkout@v4
- name: Install JuliaFormatter and format
# This will use the latest version by default but you can set the version like so:
#
diff --git a/.github/workflows/Invalidations.yml b/.github/workflows/Invalidations.yml
index 4d0004e83..28b9ce2fa 100644
--- a/.github/workflows/Invalidations.yml
+++ b/.github/workflows/Invalidations.yml
@@ -19,12 +19,12 @@ jobs:
- uses: julia-actions/setup-julia@v1
with:
version: '1'
- - uses: actions/checkout@v3
+ - uses: actions/checkout@v4
- uses: julia-actions/julia-buildpkg@v1
- uses: julia-actions/julia-invalidations@v1
id: invs_pr
- - uses: actions/checkout@v3
+ - uses: actions/checkout@v4
with:
ref: ${{ github.event.repository.default_branch }}
- uses: julia-actions/julia-buildpkg@v1
diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml
new file mode 100644
index 000000000..9246edd2a
--- /dev/null
+++ b/.github/workflows/SpellCheck.yml
@@ -0,0 +1,13 @@
+name: Spell Check
+
+on: [pull_request]
+
+jobs:
+ typos-check:
+ name: Spell Check with Typos
+ runs-on: ubuntu-latest
+ steps:
+ - name: Checkout Actions Repository
+ uses: actions/checkout@v4
+ - name: Check spelling
+ uses: crate-ci/typos@v1.18.0
\ No newline at end of file
diff --git a/.typos.toml b/.typos.toml
new file mode 100644
index 000000000..9b079dcfc
--- /dev/null
+++ b/.typos.toml
@@ -0,0 +1,4 @@
+[default.extend-words]
+Nd = "Nd"
+nin = "nin"
+coul = "coul"
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index af522a647..c4e52f67b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,25 +1,34 @@
name = "ModelingToolkitStandardLibrary"
uuid = "16a59e39-deab-5bd0-87e4-056b12336739"
authors = ["Chris Rackauckas and Julia Computing"]
-version = "2.1.1"
+version = "2.3.5"
[deps]
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
[compat]
+Aqua = "0.8"
ChainRulesCore = "1"
-DiffEqBase = "6"
+ControlSystemsBase = "1"
+DataInterpolations = "4"
+DiffEqBase = "6.103"
IfElse = "0.1"
-ModelingToolkit = "8.50"
-Symbolics = "4.9, 5"
-julia = "1.6"
+LinearAlgebra = "1.10"
+ModelingToolkit = "8.69"
+OrdinaryDiffEq = "6.33"
+SafeTestsets = "0.1"
+Symbolics = "5.2"
+Test = "1"
+julia = "1.10"
[extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -28,4 +37,4 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
[targets]
-test = ["LinearAlgebra", "OrdinaryDiffEq", "SafeTestsets", "Test", "ControlSystemsBase", "DataInterpolations"]
+test = ["Aqua", "LinearAlgebra", "OrdinaryDiffEq", "SafeTestsets", "Test", "ControlSystemsBase", "DataInterpolations"]
diff --git a/README.md b/README.md
index 52e6f84e9..94dd89676 100644
--- a/README.md
+++ b/README.md
@@ -32,7 +32,7 @@ the documentation, which contains the unreleased features.
## Libraries
-The following are the constituant libraries of the ModelingToolkit Standard Library.
+The following are the constituent libraries of the ModelingToolkit Standard Library.
- [Basic Blocks](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/blocks/)
- [Mechanical Components](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/mechanical/)
@@ -60,9 +60,9 @@ V = 1.0
@named ground = Ground()
rc_eqs = [connect(constant.output, source.V)
- connect(source.p, resistor.p)
- connect(resistor.n, capacitor.p)
- connect(capacitor.n, source.n, ground.g)]
+ 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, constant, source, ground])
diff --git a/docs/Project.toml b/docs/Project.toml
index 5a5c39a45..034ba981d 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -11,9 +11,8 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
[compat]
ControlSystemsBase = "1.1"
DifferentialEquations = "7.6"
-Documenter = "0.27"
+Documenter = "1"
IfElse = "0.1"
-ModelingToolkit = "8"
-ModelingToolkitStandardLibrary = "2"
+ModelingToolkit = "8.67"
OrdinaryDiffEq = "6.31"
Plots = "1.36"
diff --git a/docs/make.jl b/docs/make.jl
index ed1072c80..3a6e05ec8 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -19,14 +19,6 @@ include("pages.jl")
makedocs(sitename = "ModelingToolkitStandardLibrary.jl",
authors = "Julia Computing",
- clean = true,
- doctest = false,
- linkcheck = true,
- strict = [
- :linkcheck,
- :doctest,
- :example_block,
- ],
modules = [ModelingToolkit,
ModelingToolkitStandardLibrary,
ModelingToolkitStandardLibrary.Blocks,
@@ -38,6 +30,8 @@ makedocs(sitename = "ModelingToolkitStandardLibrary.jl",
ModelingToolkitStandardLibrary.Thermal,
ModelingToolkitStandardLibrary.Hydraulic,
ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible],
+ clean = true, doctest = false, linkcheck = true,
+ warnonly = [:docs_block, :missing_docs, :cross_references],
format = Documenter.HTML(assets = ["assets/favicon.ico"],
canonical = "https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/"),
pages = pages)
diff --git a/docs/pages.jl b/docs/pages.jl
index af2de24f3..35cf400cb 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -5,7 +5,7 @@ pages = [
"Custom Components" => "tutorials/custom_component.md",
"Thermal Conduction Model" => "tutorials/thermal_model.md",
"DC Motor with Speed Controller" => "tutorials/dc_motor_pi.md",
- "SampledData Component" => "tutorials/input_component.md",
+ "SampledData Component" => "tutorials/input_component.md"
],
"About Acausal Connections" => "connectors/connections.md",
"API" => [
@@ -15,6 +15,6 @@ pages = [
"Mechanical Components" => "API/mechanical.md",
"Thermal Components" => "API/thermal.md",
"Hydraulic Components" => "API/hydraulic.md",
- "Linear Analysis" => "API/linear_analysis.md",
- ],
+ "Linear Analysis" => "API/linear_analysis.md"
+ ]
]
diff --git a/docs/src/API/blocks.md b/docs/src/API/blocks.md
index 0c113fbdc..11f241fd7 100644
--- a/docs/src/API/blocks.md
+++ b/docs/src/API/blocks.md
@@ -80,6 +80,7 @@ Derivative
FirstOrder
SecondOrder
StateSpace
+TransferFunction
PI
LimPI
PID
diff --git a/docs/src/API/linear_analysis.md b/docs/src/API/linear_analysis.md
index 0028defb8..98d265371 100644
--- a/docs/src/API/linear_analysis.md
+++ b/docs/src/API/linear_analysis.md
@@ -60,7 +60,7 @@ using ModelingToolkitStandardLibrary.Blocks, ModelingToolkit
@named C = Gain(-1) # A P controller
t = ModelingToolkit.get_iv(P)
eqs = [connect(P.output, :plant_output, C.input) # Connect with an automatically created analysis point called :plant_output
- connect(C.output, :plant_input, P.input)]
+ connect(C.output, :plant_input, P.input)]
sys = ODESystem(eqs, t, systems = [P, C], name = :feedback_system)
matrices_S = get_sensitivity(sys, :plant_input)[1] # Compute the matrices of a state-space representation of the (input)sensitivity function.
diff --git a/docs/src/connectors/connections.md b/docs/src/connectors/connections.md
index 28ec45c92..ac4369977 100644
--- a/docs/src/connectors/connections.md
+++ b/docs/src/connectors/connections.md
@@ -90,16 +90,15 @@ As can be seen, this will give a 1 equation model matching our energy dissipatio
```@example connections
using ModelingToolkitStandardLibrary.Electrical, ModelingToolkit, DifferentialEquations
+using ModelingToolkit: t, D
using Plots
-@parameters t
-
@named resistor = Resistor(R = 1)
@named capacitor = Capacitor(C = 1)
@named ground = Ground()
eqs = [connect(capacitor.p, resistor.p)
- connect(resistor.n, ground.g, capacitor.n)]
+ connect(resistor.n, ground.g, capacitor.n)]
@named model = ODESystem(eqs, t; systems = [resistor, capacitor, ground])
@@ -139,7 +138,7 @@ const TV = ModelingToolkitStandardLibrary.Mechanical.Translational
@named ground = TV.Fixed()
eqs = [connect(damping.flange_a, body.flange)
- connect(ground.flange, damping.flange_b)]
+ connect(ground.flange, damping.flange_b)]
@named model = ODESystem(eqs, t; systems = [damping, body, ground])
@@ -172,7 +171,7 @@ const TP = ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition
@named ground = TP.Fixed(s_0 = 0)
eqs = [connect(damping.flange_a, body.flange)
- connect(ground.flange, damping.flange_b)]
+ connect(ground.flange, damping.flange_b)]
@named model = ODESystem(eqs, t; systems = [damping, body, ground])
@@ -267,7 +266,7 @@ Let's define a quick function to simplify and solve the 2 different systems. Not
```@example connections
function simplify_and_solve(damping, spring, body, ground)
eqs = [connect(spring.flange_a, body.flange, damping.flange_a)
- connect(spring.flange_b, damping.flange_b, ground.flange)]
+ connect(spring.flange_b, damping.flange_b, ground.flange)]
@named model = ODESystem(eqs, t; systems = [ground, body, spring, damping])
diff --git a/docs/src/index.md b/docs/src/index.md
index 90aed9b00..b0319ca60 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -21,7 +21,7 @@ Pkg.add("ModelingToolkitStandardLibrary")
## Libraries
-The following are the constituant libraries of the ModelingToolkit Standard Library.
+The following are the constituent libraries of the ModelingToolkit Standard Library.
- [Basic Blocks](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/blocks/)
- [Mechanical Components](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/mechanical/)
@@ -87,32 +87,19 @@ Pkg.status(; mode = PKGMODE_MANIFEST) # hide
```
-```@raw html
-You can also download the
-manifest file and the
-project file.
+link_manifest = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version *
+ "/assets/Manifest.toml"
+link_project = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version *
+ "/assets/Project.toml"
+Markdown.parse("""You can also download the
+[manifest]($link_manifest)
+file and the
+[project]($link_project)
+file.
+""")
```
diff --git a/docs/src/tutorials/custom_component.md b/docs/src/tutorials/custom_component.md
index 53f51c73e..49a3bce6a 100644
--- a/docs/src/tutorials/custom_component.md
+++ b/docs/src/tutorials/custom_component.md
@@ -8,6 +8,7 @@ First, we need to make some imports.
```@example components
using ModelingToolkit
+using ModelingToolkit: t_nounits as t
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Electrical: OnePort
using OrdinaryDiffEq
@@ -35,18 +36,16 @@ end NonlinearResistor;
this can almost be directly translated to the syntax of `ModelingToolkit`.
```@example components
-@parameters t
-
function NonlinearResistor(; name, Ga, Gb, Ve)
@named oneport = OnePort()
@unpack v, i = oneport
pars = @parameters Ga=Ga Gb=Gb Ve=Ve
eqs = [
i ~ ifelse(v < -Ve,
- Gb * (v + Ve) - Ga * Ve,
- ifelse(v > Ve,
- Gb * (v - Ve) + Ga * Ve,
- Ga * v)),
+ Gb * (v + Ve) - Ga * Ve,
+ ifelse(v > Ve,
+ Gb * (v - Ve) + Ga * Ve,
+ Ga * v))
]
extend(ODESystem(eqs, t, [], pars; name = name), oneport)
end
@@ -102,14 +101,14 @@ The final model can now be created with the components from the library and the
@named Gnd = Ground()
connections = [connect(L.p, G.p)
- connect(G.n, Nr.p)
- connect(Nr.n, Gnd.g)
- connect(C1.p, G.n)
- connect(L.n, Ro.p)
- connect(G.p, C2.p)
- connect(C1.n, Gnd.g)
- connect(C2.n, Gnd.g)
- connect(Ro.n, Gnd.g)]
+ connect(G.n, Nr.p)
+ connect(Nr.n, Gnd.g)
+ connect(C1.p, G.n)
+ connect(L.n, Ro.p)
+ connect(G.p, C2.p)
+ connect(C1.n, Gnd.g)
+ connect(C2.n, Gnd.g)
+ connect(Ro.n, Gnd.g)]
@named model = ODESystem(connections, t, systems = [L, Ro, G, C1, C2, Nr, Gnd])
nothing # hide
@@ -132,7 +131,7 @@ Plots.savefig("chua_phase_plane.png")
nothing # hide
Plots.plot(sol; idxs = [C1.v, C2.v, L.i],
- labels = ["C1 Voltage in V" "C1 Voltage in V" "Inductor Current in A"])
+ labels = ["C1 Voltage in V" "C2 Voltage in V" "Inductor Current in A"])
Plots.savefig("chua.png")
nothing # hide
```
diff --git a/docs/src/tutorials/dc_motor_pi.md b/docs/src/tutorials/dc_motor_pi.md
index 0dcc3af45..040e57a52 100644
--- a/docs/src/tutorials/dc_motor_pi.md
+++ b/docs/src/tutorials/dc_motor_pi.md
@@ -14,20 +14,19 @@ First, the needed packages are imported and the parameters of the model defined.
```@example dc_motor_pi
using ModelingToolkit
+using ModelingToolkit: t
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Mechanical.Rotational
using ModelingToolkitStandardLibrary.Blocks
using OrdinaryDiffEq
using Plots
-@parameters t
-
R = 0.5 # [Ohm] armature resistance
L = 4.5e-3 # [H] armature inductance
k = 0.5 # [N.m/A] motor constant
J = 0.02 # [kg.m²] inertia
f = 0.01 # [N.m.s/rad] friction factor
-tau_L_step = -0.3 # [N.m] amplitude of the load torque step
+tau_L_step = -0.3 # [N.m] amplitude of the load torque step
nothing # hide
```
@@ -43,25 +42,25 @@ The actual model can now be composed.
@named L1 = Inductor(L = L)
@named emf = EMF(k = k)
@named fixed = Fixed()
-@named load = Torque(use_support = false)
+@named load = Torque()
@named load_step = Blocks.Step(height = tau_L_step, start_time = 3)
@named inertia = Inertia(J = J)
@named friction = Damper(d = f)
@named speed_sensor = SpeedSensor()
connections = [connect(fixed.flange, emf.support, friction.flange_b)
- connect(emf.flange, friction.flange_a, inertia.flange_a)
- connect(inertia.flange_b, load.flange)
- connect(inertia.flange_b, speed_sensor.flange)
- connect(load_step.output, load.tau)
- connect(ref.output, feedback.input1)
- connect(speed_sensor.w, :y, feedback.input2)
- connect(feedback.output, pi_controller.err_input)
- connect(pi_controller.ctr_output, :u, source.V)
- connect(source.p, R1.p)
- connect(R1.n, L1.p)
- connect(L1.n, emf.p)
- connect(emf.n, source.n, ground.g)]
+ connect(emf.flange, friction.flange_a, inertia.flange_a)
+ connect(inertia.flange_b, load.flange)
+ connect(inertia.flange_b, speed_sensor.flange)
+ connect(load_step.output, load.tau)
+ connect(ref.output, feedback.input1)
+ connect(speed_sensor.w, :y, feedback.input2)
+ connect(feedback.output, pi_controller.err_input)
+ connect(pi_controller.ctr_output, :u, source.V)
+ connect(source.p, R1.p)
+ connect(R1.n, L1.p)
+ connect(L1.n, emf.p)
+ connect(emf.n, source.n, ground.g)]
@named model = ODESystem(connections, t,
systems = [
@@ -78,7 +77,7 @@ connections = [connect(fixed.flange, emf.support, friction.flange_b)
load_step,
inertia,
friction,
- speed_sensor,
+ speed_sensor
])
nothing # hide
```
diff --git a/docs/src/tutorials/input_component.md b/docs/src/tutorials/input_component.md
index cb4f50041..1a59a8dc3 100644
--- a/docs/src/tutorials/input_component.md
+++ b/docs/src/tutorials/input_component.md
@@ -14,13 +14,11 @@ The `ModelingToolkitStandardLibrary.Blocks.TimeVaryingFunction` component is eas
```julia
using ModelingToolkit
+using ModelingToolkit: t, D
using ModelingToolkitStandardLibrary.Blocks
using DataInterpolations
using OrdinaryDiffEq
-@parameters t
-D = Differential(t)
-
function System(f; name)
src = TimeVaryingFunction(f)
@@ -28,9 +26,9 @@ function System(f; name)
pars = @parameters m=10 k=1000 d=1
eqs = [f ~ src.output.u
- ddx * 10 ~ k * x + d * dx + f
- D(x) ~ dx
- D(dx) ~ ddx]
+ ddx * 10 ~ k * x + d * dx + f
+ D(x) ~ dx
+ D(dx) ~ ddx]
ODESystem(eqs, t, vars, pars; systems = [src], name)
end
@@ -47,7 +45,7 @@ prob = ODEProblem(sys, [], (0, time[end]))
sol = solve(prob, ImplicitEuler())
```
-If we want to run a new data set, this requires building a new `LinearInterpolation` and `ODESystem` followed by running `structural_simplify`, all of which takes time. Therefore, to run serveral pieces of data it's better to re-use an `ODESystem`. The next couple methods will demonstrate how to do this.
+If we want to run a new data set, this requires building a new `LinearInterpolation` and `ODESystem` followed by running `structural_simplify`, all of which takes time. Therefore, to run several pieces of data it's better to re-use an `ODESystem`. The next couple methods will demonstrate how to do this.
## Custom Component with External Data
@@ -74,9 +72,9 @@ function System(; name)
pars = @parameters m=10 k=1000 d=1
eqs = [f ~ get_sampled_data(t)
- ddx * 10 ~ k * x + d * dx + f
- D(x) ~ dx
- D(dx) ~ ddx]
+ ddx * 10 ~ k * x + d * dx + f
+ D(x) ~ dx
+ D(dx) ~ ddx]
ODESystem(eqs, t, vars, pars; name)
end
@@ -105,7 +103,7 @@ Additional code could be added to resolve this issue, for example by using a `Re
## `SampledData` Component
-To resolve the issues presented above, the `ModelingToolkitStandardLibrary.Blocks.SampledData` component can be used which allows for a resusable `ODESystem` and self contained data which ensures a solution which remains valid for it's lifetime. Now it's possible to also parallize the call to `solve()`.
+To resolve the issues presented above, the `ModelingToolkitStandardLibrary.Blocks.SampledData` component can be used which allows for a resusable `ODESystem` and self contained data which ensures a solution which remains valid for it's lifetime. Now it's possible to also parallelize the call to `solve()`.
```julia
function System(; name)
@@ -115,9 +113,9 @@ function System(; name)
pars = @parameters m=10 k=1000 d=1
eqs = [f ~ src.output.u
- ddx * 10 ~ k * x + d * dx + f
- D(x) ~ dx
- D(dx) ~ ddx]
+ ddx * 10 ~ k * x + d * dx + f
+ D(x) ~ dx
+ D(dx) ~ ddx]
ODESystem(eqs, t, vars, pars; systems = [src], name)
end
@@ -146,4 +144,4 @@ sol2 = Ref{ODESolution}()
end
```
-Note, in the above example, we can build the system with an empty `SampledData` component, only setting the expected data type: `@named src = SampledData(Float64)`. It's also possible to initialize the component with real sampled data: `@named src = SampledData(data, dt)`. Additionally note that before running an `ODEProblem` using the `SampledData` component, one must be careful about the parameter vector Type. The `SampledData` component contains a `buffer` parameter of type `Parameter`, therefore we must generate the problem using `tofloat=false`. This will initially give a parameter vector of type `Vector{Any}` with a mix of numbers and `Parameter` type. We can convert the vector to a uniform `Parameter` type by running `p = Parameter.(p)`. This will wrap all the single values in a `Parameter` which will be mathmatically equivalent to a `Number`.
+Note, in the above example, we can build the system with an empty `SampledData` component, only setting the expected data type: `@named src = SampledData(Float64)`. It's also possible to initialize the component with real sampled data: `@named src = SampledData(data, dt)`. Additionally note that before running an `ODEProblem` using the `SampledData` component, one must be careful about the parameter vector Type. The `SampledData` component contains a `buffer` parameter of type `Parameter`, therefore we must generate the problem using `tofloat=false`. This will initially give a parameter vector of type `Vector{Any}` with a mix of numbers and `Parameter` type. We can convert the vector to a uniform `Parameter` type by running `p = Parameter.(p)`. This will wrap all the single values in a `Parameter` which will be mathematically equivalent to a `Number`.
diff --git a/docs/src/tutorials/rc_circuit.md b/docs/src/tutorials/rc_circuit.md
index a33a89ca1..4004d5abd 100644
--- a/docs/src/tutorials/rc_circuit.md
+++ b/docs/src/tutorials/rc_circuit.md
@@ -22,9 +22,9 @@ V = 1.0
@named ground = Ground()
rc_eqs = [connect(constant.output, source.V)
- connect(source.p, resistor.p)
- connect(resistor.n, capacitor.p)
- connect(capacitor.n, source.n, ground.g)]
+ 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, constant, source, ground])
diff --git a/docs/src/tutorials/thermal_model.md b/docs/src/tutorials/thermal_model.md
index b83a064bb..9d104a615 100644
--- a/docs/src/tutorials/thermal_model.md
+++ b/docs/src/tutorials/thermal_model.md
@@ -8,8 +8,7 @@ from dividing the total initial energy in the system by the sum of the heat capa
```@example
using ModelingToolkitStandardLibrary.Thermal, ModelingToolkit, OrdinaryDiffEq, Plots
-
-@parameters t
+using ModelingToolkit: t
C1 = 15
C2 = 15
@@ -23,7 +22,7 @@ connections = [
connect(mass1.port, conduction.port_a),
connect(conduction.port_b, mass2.port),
connect(mass1.port, Tsensor1.port),
- connect(mass2.port, Tsensor2.port),
+ connect(mass2.port, Tsensor2.port)
]
@named model = ODESystem(connections, t,
diff --git a/src/Blocks/Blocks.jl b/src/Blocks/Blocks.jl
index dbf067e3f..5503b0b5c 100644
--- a/src/Blocks/Blocks.jl
+++ b/src/Blocks/Blocks.jl
@@ -5,10 +5,8 @@ module Blocks
using ModelingToolkit, Symbolics
import IfElse: ifelse
import ..@symcheck
-using ModelingToolkit: getdefault
-
-@parameters t
-D = Differential(t)
+using DynamicQuantities
+using ModelingToolkit: getdefault, t, D
export RealInput, RealOutput, SISO
include("utils.jl")
@@ -19,18 +17,18 @@ export Log, Log10
include("math.jl")
export Constant, TimeVaryingFunction, Sine, Cosine, ContinuousClock, Ramp, Step, ExpSine,
- Square, Triangular, Parameter, SampledData
+ Square, Triangular, Parameter, SampledData
include("sources.jl")
export Limiter, DeadZone, SlewRateLimiter
include("nonlinear.jl")
-export Integrator, Derivative, FirstOrder, SecondOrder, StateSpace
+export Integrator, Derivative, FirstOrder, SecondOrder, StateSpace, TransferFunction
export PI, LimPI, PID, LimPID
include("continuous.jl")
export AnalysisPoint, get_sensitivity, get_comp_sensitivity,
- get_looptransfer, open_loop
+ get_looptransfer, open_loop
include("analysis_points.jl")
end
diff --git a/src/Blocks/analysis_points.jl b/src/Blocks/analysis_points.jl
index a6b90866d..c5ee4b3a3 100644
--- a/src/Blocks/analysis_points.jl
+++ b/src/Blocks/analysis_points.jl
@@ -5,15 +5,25 @@ Base.@kwdef mutable struct AnalysisPoint
out = nothing
name::Symbol = :nothing
end
+Base.broadcastable(x::AnalysisPoint) = Ref(x)
+if Base.isdefined(ModelingToolkit, :isconnection)
+ ModelingToolkit.isconnection(::AnalysisPoint) = true
+end
Base.nameof(ap::AnalysisPoint) = ap.name
+function Base.hash(ap::AnalysisPoint, seed::UInt)
+ h1 = hash(ap.in, seed)
+ h2 = hash(ap.out, h1)
+ h3 = hash(ap.name, h2)
+ h3 ⊻ (0xd29cdc51aa6562d4 % UInt)
+end
function ap_var(sys)
if hasproperty(sys, :u)
# collect to turn symbolic arrays into arrays of symbols
return length(sys.u) == 1 ? sys.u : collect(sys.u)
end
- x = states(sys)
+ x = unknowns(sys)
length(x) == 1 && return x[1]
error("Could not determine the analysis-point variable in system $(nameof(sys)). To use an analysis point, apply it to a connection between two causal blocks containing connectors of type `RealInput/RealOutput` from ModelingToolkitStandardLibrary.Blocks.")
end
@@ -52,7 +62,7 @@ eqs = [connect(P.output, C.input)
connect(C.output, :plant_input, P.input)]
sys = ODESystem(eqs, t, systems = [P, C], name = :feedback_system)
-matrices_S, _ = get_sensitivity(sys, :plant_input) # Compute the matrices of a state-space representation of the (input) sensitivity funciton.
+matrices_S, _ = get_sensitivity(sys, :plant_input) # Compute the matrices of a state-space representation of the (input) sensitivity function.
matrices_T, _ = get_comp_sensitivity(sys, :plant_input)
```
@@ -255,7 +265,7 @@ end
const SymOrVec = Union{Symbol, Vector{Symbol}}
function get_sensitivity_function(sys, ap_name::SymOrVec; loop_openings = nothing,
- kwargs...)
+ kwargs...)
find = namespaced_ap_match(ap_name, loop_openings)
t = get_iv(sys)
aps = []
@@ -284,7 +294,7 @@ function get_sensitivity_function(sys, ap_name::SymOrVec; loop_openings = nothin
end
function get_comp_sensitivity_function(sys, ap_name::SymOrVec; loop_openings = nothing,
- kwargs...)
+ kwargs...)
find = namespaced_ap_match(ap_name, loop_openings)
t = get_iv(sys)
aps = []
@@ -313,7 +323,7 @@ function get_comp_sensitivity_function(sys, ap_name::SymOrVec; loop_openings = n
end
function get_looptransfer_function(sys, ap_name::SymOrVec; loop_openings = nothing,
- kwargs...)
+ kwargs...)
find = namespaced_ap_match(ap_name, loop_openings)
t = get_iv(sys)
aps = []
@@ -380,9 +390,9 @@ function open_loop(sys, ap_name::Symbol; ground_input = false, kwargs...)
end
function ModelingToolkit.linearization_function(sys::ModelingToolkit.AbstractSystem,
- input_name::SymOrVec, output_name;
- loop_openings = nothing,
- kwargs...)
+ input_name::SymOrVec, output_name;
+ loop_openings = nothing,
+ kwargs...)
t = get_iv(sys)
@variables u(t)=0 [input = true]
names = [input_name;]
@@ -417,7 +427,7 @@ function ModelingToolkit.linearization_function(sys::ModelingToolkit.AbstractSys
push!(multiplicities_y, length(yi))
append!(y, yi)
[ap_var(ap.in) .~ yi;
- ap_var(ap.out) .~ ap_var(ap.in)], yi
+ ap_var(ap.out) .~ ap_var(ap.in)], yi
else # loop opening
[ap_var(ap.out) .~ 0;], []
end
@@ -444,16 +454,23 @@ for f in [:get_sensitivity, :get_comp_sensitivity, :get_looptransfer, :open_loop
end
function ModelingToolkit.linearize(sys, input::AnalysisPoint, output::AnalysisPoint;
- kwargs...)
+ kwargs...)
ModelingToolkit.linearize(sys, nameof(input), nameof(output); kwargs...)
end
# Methods above are implemented in terms of linearization_function, the method below creates wrappers for linearize
for f in [:get_sensitivity, :get_comp_sensitivity, :get_looptransfer]
- @eval function $f(sys, ap, args...; loop_openings = nothing, kwargs...)
- lin_fun, ssys = $(Symbol(string(f) * "_function"))(sys, ap, args...; loop_openings,
+ @eval function $f(sys,
+ ap,
+ args...;
+ loop_openings = nothing,
+ op = Dict(),
+ p = DiffEqBase.NullParameters(),
+ kwargs...)
+ lin_fun, ssys = $(Symbol(string(f) * "_function"))(sys, ap, args...; op, p,
+ loop_openings,
kwargs...)
- ModelingToolkit.linearize(ssys, lin_fun; kwargs...), ssys
+ ModelingToolkit.linearize(ssys, lin_fun; op, p, kwargs...), ssys
end
end
@@ -465,7 +482,7 @@ Linearize a system between two analysis points. To get a loop-transfer function,
The output is allowed to be either an analysis-point name, or a vector of symbolic variables like the standard interface to `linearize`. The input must be an analysis-point name.
"""
function ModelingToolkit.linearize(sys, input_name::SymOrVec, output_name;
- loop_openings = nothing, kwargs...)
+ loop_openings = nothing, kwargs...)
lin_fun, ssys = linearization_function(sys, input_name, output_name; loop_openings,
kwargs...)
ModelingToolkit.linearize(ssys, lin_fun; kwargs...), ssys
@@ -513,6 +530,11 @@ get_comp_sensitivity
Compute the (linearized) loop-transfer function in analysis point `ap`, from `ap.out` to `ap.in`.
+!!! info "Negative feedback"
+
+ Feedback loops often use negative feedback, and the computed loop-transfer function will in this case have the negative feedback included. Standard analysis tools often assume a loop-transfer function without the negative gain built in, and the result of this function may thus need negation before use.
+
+
!!! danger "Experimental"
The analysis-point interface is currently experimental and at any time subject to breaking changes not respecting semantic versioning.
diff --git a/src/Blocks/continuous.jl b/src/Blocks/continuous.jl
index 852120189..1f2a3763d 100644
--- a/src/Blocks/continuous.jl
+++ b/src/Blocks/continuous.jl
@@ -1,8 +1,8 @@
"""
Integrator(;name, k = 1, x = 0.0)
-Outputs `y = ∫k*u dt`, corresponding to the transfer function `1/s`.
-Initial value of integrator state `x` can be set with `x`
+Outputs `y = ∫k*u dt`, corresponding to the transfer function ``1/s``.
+Initial value of integrator state ``x`` can be set with `x`
# Connectors:
@@ -19,7 +19,7 @@ Initial value of integrator state `x` can be set with `x`
x(t) = 0.0, [description = "State of Integrator"]
end
@parameters begin
- k = 1, [description = "Gain of Integrator"]
+ k = 1, [description = "Gain"]
end
@equations begin
D(x) ~ k * u
@@ -30,7 +30,7 @@ end
Derivative(; name, k = 1, T, x = 0.0)
Outputs an approximate derivative of the input. The transfer function of this block is
-Initial value of the state `x` can be set with `x`
+Initial value of the state ``x`` can be set with `x`
```
k k
@@ -57,11 +57,11 @@ A smaller `T` leads to a more ideal approximation of the derivative.
@mtkmodel Derivative begin
@extend u, y = siso = SISO()
@variables begin
- x(t) = 0.0, [description = "State of Derivative"]
+ x(t) = 0.0, [description = "Derivative-filter state"]
end
@parameters begin
- T = T, [description = "Time constant of Derivative"]
- k = 1, [description = "Gain of Derivative"]
+ T = T, [description = "Time constant"]
+ k = 1, [description = "Gain"]
end
begin
@symcheck T > 0 ||
@@ -77,8 +77,8 @@ end
FirstOrder(; name, k = 1.0, T, x = 0.0, lowpass = true)
A first-order filter with a single real pole in `s = -T` and gain `k`. If `lowpass=true` (default), the transfer function
-is given by `Y(s)/U(s) = `
-Initial value of the state `x` can be set with `x`
+is given by ``Y(s)/U(s) = ``
+
```
k
@@ -94,10 +94,12 @@ sT + 1 - k
sT + 1
```
+Initial value of the state `x` can be set with `x`
+
# Parameters:
- `k`: Gain
- - `T`: [s] Time constants (T>0 required)
+ - `T`: [s] Time constant (T>0 required)
# Connectors:
@@ -108,13 +110,15 @@ See also [`SecondOrder`](@ref)
"""
@mtkmodel FirstOrder begin
@extend u, y = siso = SISO()
+ @structural_parameters begin
+ lowpass = true
+ end
@variables begin
x(t) = 0.0, [description = "State of FirstOrder filter"]
end
@parameters begin
- lowpass = true
- T = T, [description = "Time constant of FirstOrder filter"]
- k = 1.0, [description = "Gain of FirstOrder"]
+ T = T, [description = "Time constant"]
+ k = 1.0, [description = "Gain"]
end
begin
@symcheck T > 0 ||
@@ -122,12 +126,12 @@ See also [`SecondOrder`](@ref)
end
@equations begin
D(x) ~ (k * u - x) / T
- getdefault(lowpass) ? y ~ x : y ~ k * u - x
+ lowpass ? y ~ x : y ~ k * u - x
end
end
"""
- SecondOrder(; name, k = 1.0, w, d, x = 0.0, xd = 0.0)
+ SecondOrder(; name, k = 1.0, w = 1.0, d = 1.0, x = 0.0, xd = 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) = `
@@ -160,9 +164,9 @@ Initial value of the state `x` can be set with `x`, and of derivative state `xd`
xd(t) = 0.0, [description = "Derivative state of SecondOrder filter"]
end
@parameters begin
- k = 1.0, [description = "Gain of SecondOrder"]
- w, [description = "Bandwidth of SecondOrder"]
- d, [description = "Relative damping of SecondOrder"]
+ k = 1.0, [description = "Gain"]
+ w = 1.0, [description = "Bandwidth (angular frequency)"]
+ d = 1.0, [description = "Relative damping"]
end
@equations begin
D(x) ~ xd
@@ -172,14 +176,19 @@ Initial value of the state `x` can be set with `x`, and of derivative state `xd`
end
"""
- PI(;name, gainPI.k = 1.0, T, int.x = 0.0)
+ PI(;name, k = 1.0, T = 1.0, int.x = 0.0)
Textbook version of a PI-controller without actuator saturation and anti-windup measure.
-The proportional gain can be set with `gainPI.k`
+The proportional gain can be set with `k`
Initial value of integrator state `x` can be set with `int.x`
-# Parameters:
+The PI controller is implemented on standard form:
+```math
+U(s) = k (1 + \\dfrac{1}{sT}) E(S)
+```
+# Parameters:
+ - `k`: Proportional gain
- `T`: [s] Integrator time constant (T>0 required)
# Connectors:
@@ -191,7 +200,8 @@ See also [`LimPI`](@ref)
"""
@mtkmodel PI begin
@parameters begin
- T, [description = "Integrator time constant"]
+ k = 1.0, [description = "Proportional gain"]
+ T = 1.0, [description = "Integrator time constant"]
end
begin
@symcheck T > 0 ||
@@ -200,7 +210,7 @@ See also [`LimPI`](@ref)
@components begin
err_input = RealInput() # control error
ctr_output = RealOutput() # control signal
- gainPI = Gain(; k = 1.0)
+ gainPI = Gain(; k)
addPI = Add()
int = Integrator(k = 1 / T, x = 0.0)
end
@@ -235,7 +245,7 @@ Text-book version of a PID-controller without actuator saturation and anti-windu
See also [`LimPID`](@ref)
"""
@component function PID(; name, k = 1, Ti = false, Td = false, Nd = 10, int__x = 0,
- der__x = 0)
+ der__x = 0)
with_I = !isequal(Ti, false)
with_D = !isequal(Td, false)
@named err_input = RealInput() # control error
@@ -272,7 +282,7 @@ See also [`LimPID`](@ref)
eqs = [
connect(err_input, addPID.input1),
connect(addPID.output, gainPID.input),
- connect(gainPID.output, ctr_output),
+ connect(gainPID.output, ctr_output)
]
if with_I
push!(eqs, connect(err_input, int.input))
@@ -294,6 +304,12 @@ end
Text-book version of a PI-controller with actuator saturation and anti-windup measure.
+The PI controller is implemented on standard form
+```math
+u(t) = sat(k (e(t) + ∫\\dfrac{1}{T}e(t) dt) )
+```
+The simplified expression above is given without the anti-windup protection.
+
# Parameters:
- `k`: Proportional gain
@@ -331,7 +347,7 @@ Text-book version of a PI-controller with actuator saturation and anti-windup me
connect(err_input, addTrack.input1),
connect(gainTrack.output, addTrack.input2),
connect(addTrack.output, int.input),
- connect(int.output, addPI.input2),
+ connect(int.output, addPI.input2)
]
ODESystem(eqs, t, [], []; name = name, systems = sys)
end
@@ -370,13 +386,13 @@ where the transfer function for the derivative includes additional filtering, se
- `ctr_output`
"""
@component 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,
- int__x = 0.0,
- der__x = 0.0)
+ Ni = Ti == 0 ? Inf : √(max(Td / Ti, 1e-6)),
+ Nd = 10,
+ u_max = Inf,
+ u_min = u_max > 0 ? -u_max : -Inf,
+ gains = false,
+ int__x = 0.0,
+ der__x = 0.0)
with_I = !isequal(Ti, false)
with_D = !isequal(Td, false)
with_AWM = Ni != Inf
@@ -443,7 +459,7 @@ where the transfer function for the derivative includes additional filtering, se
connect(addP.output, addPID.input1),
connect(addPID.output, gainPID.input),
connect(gainPID.output, limiter.input),
- connect(limiter.output, ctr_output),
+ connect(limiter.output, ctr_output)
]
if with_I
push!(eqs, connect(reference, addI.input1))
@@ -506,7 +522,7 @@ y &= h(x, u)
linearized around the operating point `x₀, u₀`, we have `y0, u0 = h(x₀, u₀), u₀`.
"""
@component function StateSpace(; A, B, C, D = nothing, x = zeros(size(A, 1)), name,
- u0 = zeros(size(B, 2)), y0 = zeros(size(C, 1)))
+ u0 = zeros(size(B, 2)), y0 = zeros(size(C, 1)))
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).")
@@ -522,18 +538,100 @@ linearized around the operating point `x₀, u₀`, we have `y0, u0 = h(x₀, u
@named input = RealInput(nin = nu)
@named output = RealOutput(nout = ny)
@variables x(t)[1:nx]=x [
- description = "State variables of StateSpace system $name",
+ description = "State variables of StateSpace system $name"
]
# 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) +
+ [D(x[i]) ~ sum(A[i, k] * x[k] for k in 1:nx) +
sum(B[i, j] * (input.u[j] - u0[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] - u0[k]) for k in 1:nu) + y0[j]
- for j in 1:ny]...,
+ for j in 1:ny]...
]
compose(ODESystem(eqs, t, vcat(x...), [], name = name), [input, output])
end
StateSpace(A, B, C, D = nothing; kwargs...) = StateSpace(; A, B, C, D, kwargs...)
+
+symbolic_eps(t) = eps(t)
+@register_symbolic symbolic_eps(t)
+
+"""
+ TransferFunction(; b, a, name)
+
+A single input, single output, linear time-invariant system provided as a transfer-function.
+```
+Y(s) = b(s) / a(s) U(s)
+```
+where `b` and `a` are vectors of coefficients of the numerator and denominator polynomials, respectively, ordered such that the coefficient of the highest power of `s` is first.
+
+The internal state realization is on controller canonical form, with state variable `x`, output variable `y` and input variable `u`. For numerical robustness, the realization used by the integrator is scaled by the last entry of the `a` parameter. The internally scaled state variable is available as `x_scaled`.
+
+To set the initial state, it's recommended to set the initial condition for `x`, and let that of `x_scaled` be computed automatically.
+
+# Parameters:
+- `b`: Numerator polynomial coefficients, e.g., `2s + 3` is specified as `[2, 3]`
+- `a`: Denominator polynomial coefficients, e.g., `s² + 2ωs + ω^2` is specified as `[1, 2ω, ω^2]`
+
+# Connectors:
+ - `input`
+ - `output`
+
+See also [`StateSpace`](@ref) which handles MIMO systems, as well as [ControlSystemsMTK.jl](https://juliacontrol.github.io/ControlSystemsMTK.jl/stable/) for an interface between [ControlSystems.jl](https://juliacontrol.github.io/ControlSystems.jl/stable/) and ModelingToolkit.jl for advanced manipulation of transfer functions and linear statespace systems. For linearization, see [`linearize`](@ref) and [Linear Analysis](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/linear_analysis/).
+"""
+@component function TransferFunction(; b = [1], a = [1, 1], name)
+ nb = length(b)
+ na = length(a)
+ nb <= na ||
+ error("Transfer function is not proper, the numerator must not be longer than the denominator")
+ nx = na - 1
+ nbb = max(0, na - nb)
+
+ @named begin
+ input = RealInput()
+ output = RealOutput()
+ end
+
+ @parameters begin
+ b[1:nb] = b,
+ [
+ description = "Numerator coefficients of transfer function (e.g., 2s + 3 is specified as [2,3])"
+ ]
+ a[1:na] = a,
+ [
+ description = "Denominator coefficients of transfer function (e.g., `s² + 2ωs + ω^2` is specified as [1, 2ω, ω^2])"
+ ]
+ bb[1:(nbb + nb)] = [zeros(nbb); b]
+ d = bb[1] / a[1], [description = "Direct feedthrough gain"]
+ end
+
+ a = collect(a)
+ @parameters a_end = ifelse(a[end] > 100 * symbolic_eps(sqrt(a' * a)), a[end], 1.0)
+
+ pars = [collect(b); a; collect(bb); d; a_end]
+ @variables begin
+ x(t)[1:nx] = zeros(nx),
+ [description = "State of transfer function on controller canonical form"]
+ x_scaled(t)[1:nx] = collect(x) * a_end, [description = "Scaled vector x"]
+ u(t), [description = "Input of transfer function"]
+ y(t), [description = "Output of transfer function"]
+ end
+
+ x = collect(x)
+ x_scaled = collect(x_scaled)
+
+ sts = [x; x_scaled; y; u]
+
+ if nx == 0
+ eqs = [y ~ d * u]
+ else
+ eqs = Equation[D(x_scaled[1]) ~ (-a[2:na]'x_scaled + a_end * u) / a[1]
+ D.(x_scaled[2:nx]) .~ x_scaled[1:(nx - 1)]
+ y ~ ((bb[2:na] - d * a[2:na])'x_scaled) / a_end + d * u
+ x .~ x_scaled ./ a_end]
+ end
+ push!(eqs, input.u ~ u)
+ push!(eqs, output.u ~ y)
+ compose(ODESystem(eqs, t, sts, pars; name = name), input, output)
+end
diff --git a/src/Blocks/math.jl b/src/Blocks/math.jl
index 86045ec52..9fb8ea3af 100644
--- a/src/Blocks/math.jl
+++ b/src/Blocks/math.jl
@@ -24,11 +24,11 @@ end
Gain.f(k; name) = Gain.f(; k, name)
"""
- MatrixGain(K::AbstractArray; name)
+ MatrixGain(; K::AbstractArray, name)
Output the product of a gain matrix with the input signal vector.
-# Parameters:
+# Structural parameters:
- `K`: Matrix gain
@@ -38,23 +38,22 @@ Output the product of a gain matrix with the input signal vector.
- `output`
"""
@mtkmodel MatrixGain begin
- @parameters begin
- K, [description = "Matrix gain"]
+ @structural_parameters begin
+ K
end
begin
- nout = size(getdefault(K), 1)
- nin = size(getdefault(K), 2)
+ nout = size(K, 1)
+ nin = size(K, 2)
end
@components begin
- input = RealInput(; nin = size(K, 2))
- output = RealOutput(; nout = size(K, 1))
+ input = RealInput(; nin = nin)
+ output = RealOutput(; nout = nout)
end
@equations begin
- [(@info i, j; output.u[i] ~ sum(getdefault(K)[i, j] * input.u[j])) for j in 1:nin
+ [output.u[i] ~ sum(K[i, j] * input.u[j] for j in 1:nin)
for i in 1:nout]...
end
end
-MatrixGain.f(K; name) = MatrixGain.f(; name, K)
"""
Sum(; input__nin::Int, name)
@@ -222,15 +221,15 @@ If the given function is not composed of simple core methods (e.g. sin, abs, ...
- `output`
"""
@mtkmodel StaticNonLinearity begin
- @parameters begin
+ @structural_parameters begin
func
end
@extend u, y = siso = SISO()
@equations begin
- y ~ first(getdefault(func))(u)
+ y ~ func(u)
end
end
-StaticNonLinearity.f(func; name) = StaticNonLinearity.f(; name = name, func = [func])
+StaticNonLinearity.f(func; name) = StaticNonLinearity(; func, name)
"""
Abs(; name)
diff --git a/src/Blocks/nonlinear.jl b/src/Blocks/nonlinear.jl
index bab526a83..8f2466f72 100644
--- a/src/Blocks/nonlinear.jl
+++ b/src/Blocks/nonlinear.jl
@@ -18,13 +18,14 @@ Limit the range of a signal.
"""
@component function Limiter(; name, y_max, y_min = y_max > 0 ? -y_max : -Inf)
@symcheck y_max ≥ y_min || throw(ArgumentError("`y_min` must be smaller than `y_max`"))
- @named siso = SISO()
+ m = (y_max + y_min) / 2
+ @named siso = SISO(u_start = m, y_start = m) # Default signals to center of saturation to minimize risk of saturation while linearizing etc.
@unpack u, y = siso
pars = @parameters y_max=y_max [description = "Maximum allowed output of Limiter $name"] y_min=y_min [
- description = "Minimum allowed output of Limiter $name",
+ description = "Minimum allowed output of Limiter $name"
]
eqs = [
- y ~ _clamp(u, y_min, y_max),
+ y ~ _clamp(u, y_min, y_max)
]
extend(ODESystem(eqs, t, [], pars; name = name), siso)
end
@@ -97,7 +98,6 @@ Initial value of state `Y` can be set with `int.y`
rising = 1.0, [description = "Maximum rising slew rate of SlewRateLimiter"]
falling = -rising, [description = "Derivative time constant of SlewRateLimiter"]
Td = 0.001, [description = "Derivative time constant"]
- y_start
end
begin
getdefault(rising) ≥ getdefault(falling) ||
@@ -105,7 +105,7 @@ Initial value of state `Y` can be set with `int.y`
getdefault(Td) > 0 ||
throw(ArgumentError("Time constant `Td` must be strictly positive"))
end
- @extend u, y = siso = SISO(y_start = y_start)
+ @extend u, y = siso = SISO(; y_start)
@equations begin
D(y) ~ max(min((u - y) / Td, rising), falling)
end
diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl
index 0b7b278a1..a1b5fda07 100644
--- a/src/Blocks/sources.jl
+++ b/src/Blocks/sources.jl
@@ -76,7 +76,7 @@ Generate constant signal.
"""
@mtkmodel Constant begin
@components begin
- output = RealOutput()
+ output = RealOutput(; unit)
end
@parameters begin
k = 0.0, [description = "Constant output value of block"]
@@ -97,17 +97,17 @@ The input variable `t` can be changed by passing a different variable as the key
- `output`
"""
@mtkmodel TimeVaryingFunction begin
- @parameters begin
+ @structural_parameters begin
f
end
@components begin
- output = RealOutput()
+ output = RealOutput(; unit)
end
@equations begin
- output.u ~ first(getdefault(f))(t)
+ output.u ~ f(t)
end
end
-TimeVaryingFunction.f(f; name) = TimeVaryingFunction.f(; f = [f], name)
+TimeVaryingFunction.f(f; name) = TimeVaryingFunction(; f, name)
"""
Sine(; name, frequency, amplitude = 1, phase = 0, offset = 0, start_time = 0,
@@ -130,13 +130,14 @@ Generate sine signal.
- `output`
"""
@component function Sine(; name,
- frequency,
- amplitude = 1,
- phase = 0,
- offset = 0,
- start_time = 0,
- smooth = false)
- @named output = RealOutput()
+ frequency,
+ amplitude = 1,
+ phase = 0,
+ offset = 0,
+ start_time = 0,
+ smooth = false,
+ output__unit = nothing)
+ @named output = RealOutput(; unit = output__unit)
pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase
equation = if smooth == false
offset + ifelse(t < start_time, 0,
@@ -147,7 +148,7 @@ Generate sine signal.
end
eqs = [
- output.u ~ equation,
+ output.u ~ equation
]
compose(ODESystem(eqs, t, [], pars; name = name), [output])
@@ -173,13 +174,14 @@ Cosine signal.
"""
@component function Cosine(; name,
- frequency,
- amplitude = 1,
- phase = 0,
- offset = 0,
- start_time = 0,
- smooth = false)
- @named output = RealOutput()
+ frequency,
+ amplitude = 1,
+ phase = 0,
+ offset = 0,
+ start_time = 0,
+ smooth = false,
+ output__unit = nothing)
+ @named output = RealOutput(; unit = output__unit)
pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase
equation = if smooth == false
offset + ifelse(t < start_time, zero(t),
@@ -189,7 +191,7 @@ Cosine signal.
smooth_cos(t, smooth, frequency, amplitude, phase, offset, start_time)
end
eqs = [
- output.u ~ equation,
+ output.u ~ equation
]
compose(ODESystem(eqs, t, [], pars; name = name), [output])
@@ -209,11 +211,11 @@ Generate current time signal.
- `output`
"""
-@component function ContinuousClock(; name, offset = 0, start_time = 0)
- @named output = RealOutput()
+@component function ContinuousClock(; name, offset = 0, start_time = 0, output__unit = nothing)
+ @named output = RealOutput(; unit = output__unit)
pars = @parameters offset=offset start_time=start_time
eqs = [
- output.u ~ offset + ifelse(t < start_time, zero(t), t - start_time),
+ output.u ~ offset + ifelse(t < start_time, zero(t), t - start_time)
]
compose(ODESystem(eqs, t, [], pars; name = name), [output])
@@ -238,12 +240,13 @@ Generate ramp signal.
- `output`
"""
@component function Ramp(; name,
- height = 1,
- duration = 1,
- offset = 0,
- start_time = 0,
- smooth = false)
- @named output = RealOutput()
+ height = 1,
+ duration = 1,
+ offset = 0,
+ start_time = 0,
+ smooth = false,
+ output__unit = nothing)
+ @named output = RealOutput(; unit = output__unit)
pars = @parameters offset=offset start_time=start_time height=height duration=duration
equation = if smooth == false
offset + ifelse(t < start_time, 0,
@@ -255,7 +258,7 @@ Generate ramp signal.
end
eqs = [
- output.u ~ equation,
+ output.u ~ equation
]
compose(ODESystem(eqs, t, [], pars; name = name), [output])
@@ -280,8 +283,9 @@ Generate smooth square signal.
- `output`
"""
@component function Square(; name, frequency = 1.0, amplitude = 1.0,
- offset = 0.0, start_time = 0.0, smooth = false)
- @named output = RealOutput()
+ offset = 0.0, start_time = 0.0, smooth = false,
+ output__unit = nothing)
+ @named output = RealOutput(; unit = output__unit)
pars = @parameters begin
frequency = frequency
amplitude = amplitude
@@ -297,7 +301,7 @@ Generate smooth square signal.
end
eqs = [
- output.u ~ equation,
+ output.u ~ equation
]
compose(ODESystem(eqs, t, [], pars; name = name), [output])
@@ -322,8 +326,8 @@ Generate step signal.
- `output`
"""
@component function Step(; name, height = 1, offset = 0, start_time = 0, duration = Inf,
- smooth = 1e-5)
- @named output = RealOutput()
+ smooth = 1e-5, output__unit = nothing)
+ @named output = RealOutput(; unit = output__unit)
duration_numeric = duration
pars = @parameters offset=offset start_time=start_time height=height duration=duration
equation = if smooth == false # use comparison in case smooth is a float
@@ -339,7 +343,7 @@ Generate step signal.
end
eqs = [
- output.u ~ equation,
+ output.u ~ equation
]
compose(ODESystem(eqs, t, [], pars; name = name), [output])
@@ -366,14 +370,15 @@ Exponentially damped sine signal.
- `output`
"""
@component function ExpSine(; name,
- frequency,
- amplitude = 1,
- damping = 0.1,
- phase = 0,
- offset = 0,
- start_time = 0,
- smooth = false)
- @named output = RealOutput()
+ frequency,
+ amplitude = 1,
+ damping = 0.1,
+ phase = 0,
+ offset = 0,
+ start_time = 0,
+ smooth = false,
+ output__unit = nothing)
+ @named output = RealOutput(; unit = output__unit)
pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase damping=damping
equation = if smooth == false
@@ -387,7 +392,7 @@ Exponentially damped sine signal.
end
eqs = [
- output.u ~ equation,
+ output.u ~ equation
]
compose(ODESystem(eqs, t, [], pars; name = name), [output])
@@ -413,8 +418,9 @@ Generate smooth triangular signal for frequencies less than or equal to 25 Hz
- `output`
"""
@component function Triangular(; name, amplitude = 1.0, frequency = 1.0,
- offset = 0.0, start_time = 0.0, smooth = false)
- @named output = RealOutput()
+ offset = 0.0, start_time = 0.0, smooth = false,
+ output__unit = nothing)
+ @named output = RealOutput(; unit = output__unit)
pars = @parameters begin
amplitude = amplitude
frequency = frequency
@@ -430,7 +436,7 @@ Generate smooth triangular signal for frequencies less than or equal to 25 Hz
end
eqs = [
- output.u ~ equation,
+ output.u ~ equation
]
compose(ODESystem(eqs, t, [], pars; name = name), [output])
@@ -442,18 +448,9 @@ end
# - SawTooth Generate saw tooth signal
# - Trapezoid Generate trapezoidal signal of type Real
-function linear_interpolation(x1::T, x2::T, t1::T, t2::T, t) where {T <: Real}
- if t1 != t2
- slope = (x2 - x1) / (t2 - t1)
- intercept = x1 - slope * t1
- return slope * t + intercept
- else
- @assert x1==x2 "x1 ($x1) and x2 ($x2) should be equal if t1 == t2"
- return x2
- end
-end
+# SampledData Parameter struct ----------------
struct Parameter{T <: Real}
data::Vector{T}
@@ -531,12 +528,55 @@ function Base.show(io::IO, m::MIME"text/plain", p::Parameter)
end
end
-function get_sampled_data(t, memory::Parameter{T}) where {T}
+get_sample_time(memory::Parameter) = memory.ref
+Symbolics.@register_symbolic get_sample_time(memory)
+
+Base.convert(::Type{T}, x::Parameter{T}) where {T <: Real} = x.ref
+function Base.convert(::Type{<:Parameter{T}}, x::Number) where {T <: Real}
+ Parameter{T}(T[], x, true)
+end
+
+
+
+# SampledData utilities ----------------
+
+function linear_interpolation(x1::Real, x2::Real, t1::Real, t2::Real, t)
+ if t1 != t2
+ slope = (x2 - x1) / (t2 - t1)
+ intercept = x1 - slope * t1
+
+ return slope * t + intercept
+ else
+ @assert x1==x2 "x1 ($x1) and x2 ($x2) should be equal if t1 == t2"
+
+ return x2
+ end
+end
+
+function first_order_backwards_difference(t, memory)
+ Δt = get_sample_time(memory)
+ x1 = get_sampled_data(t, memory)
+ x0 = get_sampled_data(t - Δt, memory)
+
+ return (x1 - x0) / Δt
+end
+
+function first_order_backwards_difference(t, buffer, Δt, circular_buffer)
+ x1 = get_sampled_data(t , buffer, Δt, circular_buffer)
+ x0 = get_sampled_data(t - Δt, buffer, Δt, circular_buffer)
+
+ return (x1 - x0) / Δt
+end
+
+function get_sampled_data(t,
+ buffer::Vector{T},
+ dt::T,
+ circular_buffer = true) where {T <: Real}
if t < 0
t = zero(t)
end
- if isempty(memory.data)
+ if isempty(buffer)
if T <: AbstractFloat
return T(NaN)
else
@@ -544,18 +584,18 @@ function get_sampled_data(t, memory::Parameter{T}) where {T}
end
end
- i1 = floor(Int, t / memory.ref) + 1 #expensive
+ i1 = floor(Int, t / dt) + 1 #expensive
i2 = i1 + 1
- t1 = (i1 - 1) * memory.ref
- x1 = @inbounds memory.data[i1]
+ t1 = (i1 - 1) * dt
+ x1 = @inbounds buffer[i1]
if t == t1
return x1
else
- n = length(memory.data)
+ n = length(buffer)
- if memory.circular_buffer
+ if circular_buffer
i1 = (i1 - 1) % n + 1
i2 = (i2 - 1) % n + 1
else
@@ -565,122 +605,129 @@ function get_sampled_data(t, memory::Parameter{T}) where {T}
end
end
- t2 = (i2 - 1) * memory.ref
- x2 = @inbounds memory.data[i2]
+ t2 = (i2 - 1) * dt
+ x2 = @inbounds buffer[i2]
return linear_interpolation(x1, x2, t1, t2, t)
end
end
+get_sampled_data(t, buffer) = get_sampled_data(t, buffer.data, buffer.ref, buffer.circular_buffer)
+Symbolics.@register_symbolic get_sampled_data(t, buffer)
+Symbolics.@register_symbolic get_sampled_data(t, buffer, dt, circular_buffer) false
-get_sample_time(memory::Parameter) = memory.ref
-Symbolics.@register_symbolic get_sample_time(memory)
-
-Symbolics.@register_symbolic get_sampled_data(t, memory)
-
-function first_order_backwards_difference(t, memory)
- Δt = get_sample_time(memory)
- x1 = get_sampled_data(t, memory)
- x0 = get_sampled_data(t - Δt, memory)
-
- return (x1 - x0) / Δt
+function Symbolics.derivative(::typeof(get_sampled_data), args::NTuple{2, Any}, ::Val{1})
+ t = @inbounds args[1]
+ buffer = @inbounds args[2]
+ first_order_backwards_difference(t, buffer)
+end
+function ChainRulesCore.frule((_, ẋ, _), ::typeof(get_sampled_data), t, buffer)
+ first_order_backwards_difference(t, buffer) * ẋ
end
-function Symbolics.derivative(::typeof(get_sampled_data), args::NTuple{2, Any}, ::Val{1})
+function Symbolics.derivative(::typeof(get_sampled_data), args::NTuple{4, Any}, ::Val{1})
t = @inbounds args[1]
- memory = @inbounds args[2]
- first_order_backwards_difference(t, memory)
+ buffer = @inbounds args[2]
+ sample_time = @inbounds args[3]
+ circular_buffer = @inbounds args[4]
+ first_order_backwards_difference(t, buffer, sample_time, circular_buffer)
end
-function ChainRulesCore.frule((_, ẋ, _), ::typeof(get_sampled_data), t, memory)
- first_order_backwards_difference(t, memory) * ẋ
+function ChainRulesCore.frule((_, ẋ, _),
+ ::typeof(get_sampled_data),
+ t,
+ buffer,
+ sample_time,
+ circular_buffer)
+ first_order_backwards_difference(t, buffer, sample_time, circular_buffer) * ẋ
+end
+
+
+
+# SampledData component ----------------
+
+module SampledDataType
+@enum Option vector_based struct_based
end
"""
- SampledData(; name, buffer)
+ SampledData(; name, buffer, sample_time, circular_buffer=true)
data input component.
# Parameters:
- - `buffer`: a `Parameter` type which holds the data and sample time
+ - `buffer::Vector{Real}`: holds the data sampled at `sample_time`
+ - `sample_time::Real`
+ - `circular_buffer::Bool = true`: how to handle `t > length(buffer)*sample_time`. If true data is considered circular, otherwise last data point is held.
# Connectors:
- `output`
"""
-@component function SampledData(; name, buffer)
+@component function SampledData(::Val{SampledDataType.vector_based};
+ name,
+ buffer,
+ sample_time,
+ circular_buffer = true)
pars = @parameters begin
- buffer = buffer
+ buffer = buffer #::Vector{Real}
+ sample_time = sample_time #::Real
+ circular_buffer = circular_buffer #::Bool
end
vars = []
systems = @named begin
output = RealOutput()
end
eqs = [
- output.u ~ get_sampled_data(t, buffer),
+ output.u ~ get_sampled_data(t, buffer, sample_time, circular_buffer)
]
return ODESystem(eqs, t, vars, pars; name, systems,
- defaults = [output.u => get_sampled_data(0.0, buffer)])
+ defaults = [
+ output.u => get_sampled_data(0.0, buffer, sample_time, circular_buffer)
+ ])
end
-@deprecate Input SampledData
-function SampledData(T::Type, circular_buffer = true; name)
- SampledData(T[], zero(T), circular_buffer; name)
-end
-function SampledData(dt::T, circular_buffer = true) where {T <: Real}
- SampledData(T[], dt, circular_buffer; name)
-end
-function SampledData(data::Vector{T}, dt::T, circular_buffer = true; name) where {T <: Real}
- SampledData(; name, buffer = Parameter(data, dt, circular_buffer))
-end
+"""
+ SampledData(; name, buffer)
-Base.convert(::Type{T}, x::Parameter{T}) where {T <: Real} = x.ref
-function Base.convert(::Type{<:Parameter{T}}, x::Number) where {T <: Real}
- Parameter{T}(T[], x, true)
-end
+data input component.
-# Beta Code for potential AE Hack ----------------------
-function set_sampled_data!(memory::Parameter{T}, t, x, Δt::Parameter{T}) where {T}
- if t < 0
- t = zero(t)
- end
+# Parameters:
+ - `buffer`: a `Parameter` type which holds the data and sample time
- if t == zero(t)
- empty!(memory.data)
+# Connectors:
+ - `output`
+"""
+@component function SampledData(::Val{SampledDataType.struct_based}; name, buffer, unit = nothing)
+ pars = @parameters begin
+ buffer = buffer #::Parameter
end
-
- n = length(memory.data)
- i = round(Int, t / Δt) + 1 #expensive
- if i == n + 1
- push!(memory.data, DiffEqBase.value(x))
- elseif i <= n
- @inbounds memory.data[i] = DiffEqBase.value(x)
- else
- error("Memory buffer skipped a step: n=$n, i=$i")
+ vars = []
+ systems = @named begin
+ output = RealOutput(; unit)
end
+ eqs = [
+ output.u ~ get_sampled_data(t, buffer)
+ ]
+ return ODESystem(eqs, t, vars, pars; name, systems,
+ defaults = [output.u => get_sampled_data(0.0, buffer)])
+end
- # memory.ref = Δt
+SampledData(x::SampledDataType.Option; kwargs...) = SampledData(Val(x); kwargs...)
- return x
-end
-Symbolics.@register_symbolic set_sampled_data!(memory, t, x, Δt)
+# struct_based
+SampledData(T::Type, circular_buffer = true; name) = SampledData(SampledDataType.struct_based; name, buffer = Parameter(T[], zero(T), circular_buffer))
-function Symbolics.derivative(::typeof(set_sampled_data!), args::NTuple{4, Any}, ::Val{2})
- memory = @inbounds args[1]
- t = @inbounds args[2]
- x = @inbounds args[3]
- Δt = @inbounds args[4]
- first_order_backwards_difference(t, x, Δt, memory)
+# vector_based
+function SampledData(sample_time::T, circular_buffer = true; name) where {T <: Real}
+ SampledData(SampledDataType.vector_based;
+ name,
+ buffer = T[],
+ sample_time,
+ circular_buffer)
end
-Symbolics.derivative(::typeof(set_sampled_data!), args::NTuple{4, Any}, ::Val{3}) = 1 #set_sampled_data returns x, therefore d/dx (x) = 1
-function ChainRulesCore.frule((_, _, ṫ, ẋ, _),
- ::typeof(set_sampled_data!),
- memory,
- t,
- x,
- Δt)
- first_order_backwards_difference(t, x, Δt, memory) * ṫ + ẋ
+function SampledData(buffer::Vector{<:Real},
+ sample_time::Real,
+ circular_buffer = true;
+ name)
+ SampledData(SampledDataType.vector_based; name, buffer, sample_time, circular_buffer)
end
-
-function first_order_backwards_difference(t, x, Δt, memory)
- x1 = set_sampled_data!(memory, t, x, Δt)
- x0 = get_sampled_data(t - Δt, memory)
-
- return (x1 - x0) / Δt
+function SampledData(; name, buffer, sample_time, circular_buffer)
+ SampledData(SampledDataType.vector_based; name, buffer, sample_time, circular_buffer)
end
diff --git a/src/Blocks/utils.jl b/src/Blocks/utils.jl
index 59db3655b..dc044a4af 100644
--- a/src/Blocks/utils.jl
+++ b/src/Blocks/utils.jl
@@ -1,14 +1,30 @@
-@connector function RealInput(; name, nin = 1, u_start = nin > 1 ? zeros(nin) : 0.0)
+@connector function RealInput(; name, nin = 1, u_start = nin > 1 ? zeros(nin) : 0.0, unit = nothing)
if nin == 1
- @variables u(t)=u_start [
- input = true,
- description = "Inner variable in RealInput $name",
- ]
+ if unit !== nothing
+ @variables u(t)=u_start [
+ input = true,
+ description = "Inner variable in RealInput $name",
+ unit = unit
+ ]
+ else
+ @variables u(t)=u_start [
+ input = true,
+ description = "Inner variable in RealInput $name",
+ ]
+ end
else
- @variables u(t)[1:nin]=u_start [
- input = true,
- description = "Inner variable in RealInput $name",
- ]
+ if unit !== nothing
+ @variables u(t)[1:nin]=u_start [
+ input = true,
+ description = "Inner variable in RealInput $name",
+ unit = unit
+ ]
+ else
+ @variables u(t)[1:nin]=u_start [
+ input = true,
+ description = "Inner variable in RealInput $name",
+ ]
+ end
u = collect(u)
end
ODESystem(Equation[], t, [u...], []; name = name)
@@ -26,18 +42,34 @@ Connector with one input signal of type Real.
- `u`: Value of the connector; if nin=1 this is a scalar
""" RealInput
-@connector function RealOutput(; name, nout = 1, u_start = nout > 1 ? zeros(nout) : 0.0)
+@connector function RealOutput(; name, nout = 1, u_start = nout > 1 ? zeros(nout) : 0.0, unit = nothing)
if nout == 1
- @variables u(t)=u_start [
- output = true,
- description = "Inner variable in RealOutput $name",
- ]
+ if unit !== nothing
+ @variables u(t)=u_start [
+ output = true,
+ description = "Inner variable in RealOutput $name",
+ unit = unit
+ ]
+ else
+ @variables u(t)=u_start [
+ output = true,
+ description = "Inner variable in RealOutput $name",
+ ]
+ end
else
- @variables u(t)[1:nout]=u_start [
- output = true,
- description = "Inner variable in RealOutput $name",
- ]
- u = collect(u)
+ if unit !== nothing
+ @variables u(t)[1:nout]=u_start [
+ output = true,
+ description = "Inner variable in RealOutput $name",
+ unit = unit
+ ]
+ else
+ @variables u(t)[1:nout]=u_start [
+ output = true,
+ description = "Inner variable in RealOutput $name",
+ ]
+ end
+ u = collect(u)
end
ODESystem(Equation[], t, [u...], []; name = name)
end
@@ -74,8 +106,8 @@ Single input single output (SISO) continuous system block.
y(t) = y_start, [description = "Output of SISO system"]
end
@components begin
- input = RealInput(u_start = 0.0)
- output = RealOutput(u_start = 0.0)
+ input = RealInput(u_start = u_start)
+ output = RealOutput(u_start = y_start)
end
@equations begin
u ~ input.u
@@ -96,14 +128,14 @@ Base class for a multiple input multiple output (MIMO) continuous system block.
- `y_start`: Initial value for the output
"""
@component function MIMO(; name, nin = 1, nout = 1, u_start = zeros(nin),
- y_start = zeros(nout))
+ y_start = zeros(nout))
@named input = RealInput(nin = nin, u_start = u_start)
@named output = RealOutput(nout = nout, u_start = y_start)
@variables(u(t)[1:nin]=u_start, [description = "Input of MIMO system $name"],
y(t)[1:nout]=y_start, [description = "Output of MIMO system $name"],)
eqs = [
[u[i] ~ input.u[i] for i in 1:nin]...,
- [y[i] ~ output.u[i] for i in 1:nout]...,
+ [y[i] ~ output.u[i] for i in 1:nout]...
]
return ODESystem(eqs, t, vcat(u..., y...), []; name = name, systems = [input, output])
end
diff --git a/src/Electrical/Analog/ideal_components.jl b/src/Electrical/Analog/ideal_components.jl
index 82a787e90..59b98cdda 100644
--- a/src/Electrical/Analog/ideal_components.jl
+++ b/src/Electrical/Analog/ideal_components.jl
@@ -38,7 +38,7 @@ See [OnePort](@ref)
@mtkmodel Resistor begin
@extend v, i = oneport = OnePort()
@parameters begin
- R, [description = "Resistance"]
+ R, [description = "Resistance", unit = u"Ω"]
end
@equations begin
v ~ i * R
@@ -66,7 +66,7 @@ See [OnePort](@ref)
@mtkmodel Conductor begin
@extend v, i = oneport = OnePort()
@parameters begin
- G, [description = "Conductance"]
+ G, [description = "Conductance", unit = u"S"]
end
@equations begin
i ~ v * G
@@ -93,13 +93,10 @@ See [OnePort](@ref)
- `C`: [`F`] Capacitance
"""
@mtkmodel Capacitor begin
+ @extend v, i = oneport = OnePort(; v = 0.0)
@parameters begin
- C, [description = "Capacitance"]
+ C, [description = "Capacitance",unit = u"F"]
end
- @variables begin
- v
- end
- @extend v, i = oneport = OnePort(; v = v)
@equations begin
D(v) ~ i / C
end
@@ -125,13 +122,10 @@ See [OnePort](@ref)
- `L`: [`H`] Inductance
"""
@mtkmodel Inductor begin
+ @extend v, i = oneport = OnePort(; i = 0.0)
@parameters begin
- L, [description = "Inductance"]
+ L, [description = "Inductance", unit = u"H"]
end
- @variables begin
- i
- end
- @extend v, i = oneport = OnePort(; i = i)
@equations begin
D(i) ~ 1 / L * v
end
@@ -211,9 +205,9 @@ Temperature dependent electrical resistor
heat_port = HeatPort()
end
@parameters begin
- R_ref = 1.0, [description = "Reference resistance"]
- T_ref = 300.15, [description = "Reference temperature"]
- alpha = 0, [description = "Temperature coefficient of resistance"]
+ R_ref = 1.0, [description = "Reference resistance", unit = u"Ω"]
+ T_ref = 300.15, [description = "Reference temperature", unit = u"K"]
+ alpha = 0, [description = "Temperature coefficient of resistance", unit = u"1/K"]
end
@variables begin
R(t) = R_ref
@@ -248,25 +242,24 @@ Electromotoric force (electric/mechanic transformer)
- `k`: [`N⋅m/A`] Transformation coefficient
"""
@mtkmodel EMF begin
+ @extend OnePort()
@components begin
- p = Pin()
- n = Pin()
flange = Flange()
support = Support()
end
@parameters begin
- k, [description = "Transformation coefficient"]
+ k, [description = "Transformation coefficient", unit = u"N*m/A"]
end
@variables begin
- v(t) = 0.0
- i(t) = 0.0
- phi(t) = 0.0
- w(t) = 0.0
+ phi(t) = 0.0, [description = "Rotation angle", unit = u"rad"]
+ w(t) = 0.0, [description = "Angular velocity", unit = u"rad/s"]
+ end
+ @extend v, i = oneport = OnePort()
+ @components begin
+ flange = Flange()
+ support = Support()
end
@equations begin
- v ~ p.v - n.v
- 0 ~ p.i + n.i
- i ~ p.i
phi ~ flange.phi - support.phi
D(phi) ~ w
k * w ~ v
diff --git a/src/Electrical/Analog/sensors.jl b/src/Electrical/Analog/sensors.jl
index 5c6e7f16d..94af63487 100644
--- a/src/Electrical/Analog/sensors.jl
+++ b/src/Electrical/Analog/sensors.jl
@@ -19,7 +19,7 @@ an ideal ammeter.
n = Pin()
end
@variables begin
- i(t)
+ i(t), [description = "Measured current", unit = u"A"]
end
@equations begin
p.v ~ n.v
@@ -46,7 +46,7 @@ Creates a circuit component which measures the potential at a pin.
p = Pin()
end
@variables begin
- phi(t)
+ phi(t), [description = "Measured potential", unit = u"V"]
end
@equations begin
p.i ~ 0
@@ -74,7 +74,7 @@ Creates a circuit component that measures the voltage across it. Analogous to an
n = Pin()
end
@variables begin
- v(t)
+ v(t), [description = "Voltage difference from positive to negative pin", unit = u"V"]
end
@equations begin
p.i ~ 0
@@ -112,7 +112,7 @@ consumed by a circuit.
current_sensor = CurrentSensor()
end
@variables begin
- power(t)
+ power(t), [description = "Power being consumed", unit = u"W"]
end
@equations begin
connect(voltage_sensor.p, pv)
@@ -150,8 +150,8 @@ Combines a [`VoltageSensor`](@ref) and a [`CurrentSensor`](@ref).
current_sensor = CurrentSensor()
end
@variables begin
- i(t) = 1.0
- v(t) = 1.0
+ i(t) = 1.0, [description = "Current", unit = u"A"]
+ v(t) = 1.0, [description = "Voltage", unit = u"V"]
end
@equations begin
connect(voltage_sensor.p, pv)
diff --git a/src/Electrical/Analog/sources.jl b/src/Electrical/Analog/sources.jl
index 807fa91f1..4d3f20b3e 100644
--- a/src/Electrical/Analog/sources.jl
+++ b/src/Electrical/Analog/sources.jl
@@ -16,7 +16,7 @@ See [OnePort](@ref)
@mtkmodel Voltage begin
@extend v, i = oneport = OnePort()
@components begin
- V = RealInput()
+ V = RealInput(unit = u"V")
end
@equations begin
v ~ V.u
@@ -41,7 +41,7 @@ See [OnePort](@ref)
@mtkmodel Current begin
@extend v, i = oneport = OnePort()
@components begin
- I = RealInput()
+ I = RealInput(unit = u"A")
end
@equations begin
i ~ I.u
diff --git a/src/Electrical/Digital/components.jl b/src/Electrical/Digital/components.jl
index ce9f31211..c3252f57d 100644
--- a/src/Electrical/Digital/components.jl
+++ b/src/Electrical/Digital/components.jl
@@ -30,9 +30,9 @@ Takes two bits as input, and outputs the sum and the carry
@variables sum(t), carry(t)
eqs = [y1.val ~ _xor(x1.val, x2.val)
- y2.val ~ _and(x1.val, x2.val)
- sum ~ y1.val
- carry ~ y2.val]
+ y2.val ~ _and(x1.val, x2.val)
+ sum ~ y1.val
+ carry ~ y2.val]
ODESystem(eqs, t, [sum, carry], [], systems = [x1, x2, y1, y2], name = name)
end
@@ -68,9 +68,9 @@ Takes three bits as input, and outputs the sum and the carry
@variables sum(t), carry(t)
eqs = [y1.val ~ _xor(x1.val, x2.val, x3.val)
- y2.val ~ _or(_and(x3.val, _xor(x1.val, x2.val)), _and(x1.val, x2.val))
- sum ~ y1.val
- carry ~ y2.val]
+ y2.val ~ _or(_and(x3.val, _xor(x1.val, x2.val)), _and(x1.val, x2.val))
+ sum ~ y1.val
+ carry ~ y2.val]
ODESystem(eqs, t, [sum, carry], [], systems = [x1, x2, x3, y1, y2], name = name)
end
diff --git a/src/Electrical/Digital/gates.jl b/src/Electrical/Digital/gates.jl
index e7cdd73f5..0940e4cfc 100644
--- a/src/Electrical/Digital/gates.jl
+++ b/src/Electrical/Digital/gates.jl
@@ -17,7 +17,7 @@ function Not(; name)
@named y = DigitalPin()
eqs = [x.i ~ y.i
- y.val ~ _not(x.val)]
+ y.val ~ _not(x.val)]
ODESystem(eqs, t, [], [], systems = [x, y], name = name)
end
@@ -43,7 +43,7 @@ function And(; name, N = 2)
vals = [k.val for k in x]
eqs = [y.val ~ _and(vals...)
- y.i ~ sum(k -> k.i, x)]
+ y.i ~ sum(k -> k.i, x)]
ODESystem(eqs, t, [], [], systems = [x..., y], name = name)
end
@@ -69,7 +69,7 @@ function Nand(; name, N = 2)
vlist = [k.val for k in x]
eqs = [y.val ~ _not(_and(vlist...))
- y.i ~ sum(k -> k.i, x)]
+ y.i ~ sum(k -> k.i, x)]
ODESystem(eqs, t, [], [], systems = [x..., y], name = name)
end
@@ -95,7 +95,7 @@ function Or(; name, N = 2)
vals = [k.val for k in x]
eqs = [y.val ~ _or(vals...)
- y.i ~ sum(k -> k.i, x)]
+ y.i ~ sum(k -> k.i, x)]
ODESystem(eqs, t, [], [], systems = [x..., y], name = name)
end
@@ -121,7 +121,7 @@ function Nor(; name, N = 2)
vlist = [k.val for k in x]
eqs = [y.val ~ _not(_or(vlist...))
- y.i ~ sum(k -> k.i, x)]
+ y.i ~ sum(k -> k.i, x)]
ODESystem(eqs, t, [], [], systems = [x..., y], name = name)
end
@@ -147,7 +147,7 @@ function Xor(; name, N = 2)
vals = [k.val for k in x]
eqs = [y.val ~ _xor(vals...)
- y.i ~ sum(k -> k.i, x)]
+ y.i ~ sum(k -> k.i, x)]
ODESystem(eqs, t, [], [], systems = [x..., y], name = name)
end
@@ -173,6 +173,6 @@ function Xnor(; name, N = 2)
vlist = [k.val for k in x]
eqs = [y.val ~ _not(_xor(vlist...))
- y.i ~ sum(k -> k.i, x)]
+ y.i ~ sum(k -> k.i, x)]
ODESystem(eqs, t, [], [], systems = [x..., y], name = name)
end
diff --git a/src/Electrical/Digital/logic.jl b/src/Electrical/Digital/logic.jl
index 675a0d7ce..ed6f52ffa 100644
--- a/src/Electrical/Digital/logic.jl
+++ b/src/Electrical/Digital/logic.jl
@@ -1,10 +1,10 @@
-@enum Logic Uninitialized=1 ForcingUnknown ForcingZero ForcingOne HighImpedence WeakUnknown WeakZero WeakOne DontCare
+@enum Logic Uninitialized=1 ForcingUnknown ForcingZero ForcingOne HighImpedance WeakUnknown WeakZero WeakOne DontCare
const U = Uninitialized
const X = ForcingUnknown
const F0 = ForcingZero
const F1 = ForcingOne
-const Z = HighImpedence
+const Z = HighImpedance
const W = WeakUnknown
const L = WeakZero
const H = WeakOne
diff --git a/src/Electrical/Digital/logic_vectors.jl b/src/Electrical/Digital/logic_vectors.jl
index 8f048249a..590621cd1 100644
--- a/src/Electrical/Digital/logic_vectors.jl
+++ b/src/Electrical/Digital/logic_vectors.jl
@@ -24,7 +24,7 @@ axes(l::LogicVector) = axes(l.logic)
getindex(s::LogicVector, i::Int) = getindex(s.logic, i)
function Base.getindex(s::LogicVector, i1::Int, i2::Int,
- I::Int...)
+ I::Int...)
getindex(s.logic, i1, i2, I...)
end
diff --git a/src/Electrical/Digital/sources.jl b/src/Electrical/Digital/sources.jl
index 17537c088..b2a41a8e8 100644
--- a/src/Electrical/Digital/sources.jl
+++ b/src/Electrical/Digital/sources.jl
@@ -19,7 +19,7 @@ function PulseDiff(; name, Val = 1, dt = 0.1)
D = ModelingToolkit.Difference(t; dt = dt)
eqs = [D(val) ~ Val
- val ~ d.val]
+ val ~ d.val]
ODESystem(eqs, t, [val], [], systems = [d], defaults = Dict(Val => 0), name = name)
end
@@ -40,7 +40,7 @@ function Set(; name)
@named d = DigitalPin()
eqs = [
- d.val ~ 1,
+ d.val ~ 1
]
ODESystem(eqs, t, [], [], systems = [d], name = name)
end
@@ -61,7 +61,7 @@ function Reset(; name)
@named d = DigitalPin()
eqs = [
- d.val ~ 0,
+ d.val ~ 0
]
ODESystem(eqs, t, [], [], systems = [d], name = name)
end
@@ -82,7 +82,7 @@ function Pulse(; name, duty_cycle = 0.5, T = 1.0)
@named d = DigitalPin()
eqs = [
- d.val ~ IfElse.ifelse(t % T > duty_cycle * T, 1, 0),
+ d.val ~ IfElse.ifelse(t % T > duty_cycle * T, 1, 0)
]
ODESystem(eqs, t, [], [], systems = [d], name = name)
end
diff --git a/src/Electrical/Digital/tables.jl b/src/Electrical/Digital/tables.jl
index b2d7b8983..beb171730 100644
--- a/src/Electrical/Digital/tables.jl
+++ b/src/Electrical/Digital/tables.jl
@@ -17,7 +17,7 @@ function Base.getindex(s::LogicTable, i1::Logic, i2::Logic)
getindex(s.logic, get_logic_level(i1), get_logic_level(i2))
end
function Base.getindex(s::LogicTable, i1::Logic, i2::Logic,
- I::Logic...)
+ I::Logic...)
getindex(s.logic, get_logic_level(i1), get_logic_level(i2), get_logic_level(I...)...)
end
@@ -37,16 +37,16 @@ get_logic_level(l::LogicTable) = Int.(l.logic)
# AND gate
const AndTable = LogicTable([
-# U X F0 F1 Z W L H DC
- U U F0 U U U F0 U U # U
- U X F0 X X X F0 X X # X
- F0 F0 F0 F0 F0 F0 F0 F0 F0 # F0
- U X F0 F1 X X F0 F1 X # F1
- U X F0 X X X F0 X X # Z
- U X F0 X X X F0 X X # W
- F0 F0 F0 F0 F0 F0 F0 F0 F0 # L
- U X F0 F1 X X F0 F1 X # H
- U X F0 X X X F0 X X]) # DC
+ # U X F0 F1 Z W L H DC
+ U U F0 U U U F0 U U # U
+ U X F0 X X X F0 X X # X
+ F0 F0 F0 F0 F0 F0 F0 F0 F0 # F0
+ U X F0 F1 X X F0 F1 X # F1
+ U X F0 X X X F0 X X # Z
+ U X F0 X X X F0 X X # W
+ F0 F0 F0 F0 F0 F0 F0 F0 F0 # L
+ U X F0 F1 X X F0 F1 X # H
+ U X F0 X X X F0 X X]) # DC
function _and2(a::Logic, b::Logic)
AndTable[a, b]
@@ -75,16 +75,16 @@ _not(x::Number) = _not(convert(Logic, x))
# OR gate
const OrTable = LogicTable([
-# U X F0 F1 Z W L H DC
- U U U F1 U U U F1 U # U
- U X X F1 X X X F1 X # X
- U X F0 F1 X X F0 F1 X # F0
- F1 F1 F1 F1 F1 F1 F1 F1 F1 # F1
- U X X F1 X X X F1 X # Z
- U X X F1 X X X F1 X # W
- U X F0 F1 X X F0 F1 X # L
- F1 F1 F1 F1 F1 F1 F1 F1 F1 # H
- U X X F1 X X X F1 X]) # DC
+ # U X F0 F1 Z W L H DC
+ U U U F1 U U U F1 U # U
+ U X X F1 X X X F1 X # X
+ U X F0 F1 X X F0 F1 X # F0
+ F1 F1 F1 F1 F1 F1 F1 F1 F1 # F1
+ U X X F1 X X X F1 X # Z
+ U X X F1 X X X F1 X # W
+ U X F0 F1 X X F0 F1 X # L
+ F1 F1 F1 F1 F1 F1 F1 F1 F1 # H
+ U X X F1 X X X F1 X]) # DC
function _or2(a::Logic, b::Logic)
OrTable[a, b]
@@ -105,16 +105,16 @@ end
# XOR gate
const XorTable = LogicTable([
-# U X F0 F1 Z W L H DC
- U U U U U U U U U # U
- U X X X X X X X X # X
- U X F0 F1 X X F0 F1 X # F0
- U X F1 F0 X X F1 F0 X # F1
- U X X X X X X X X # Z
- U X X X X X X X X # W
- U X F0 F1 X X F0 F1 X # L
- U X F1 F0 X X F1 F0 X # H
- U X X X X X X X X]) # DC
+ # U X F0 F1 Z W L H DC
+ U U U U U U U U U # U
+ U X X X X X X X X # X
+ U X F0 F1 X X F0 F1 X # F0
+ U X F1 F0 X X F1 F0 X # F1
+ U X X X X X X X X # Z
+ U X X X X X X X X # W
+ U X F0 F1 X X F0 F1 X # L
+ U X F1 F0 X X F1 F0 X # H
+ U X X X X X X X X]) # DC
function _xor2(a::Logic, b::Logic)
XorTable[a, b]
diff --git a/src/Electrical/Electrical.jl b/src/Electrical/Electrical.jl
index 687f92b20..847ef2470 100644
--- a/src/Electrical/Electrical.jl
+++ b/src/Electrical/Electrical.jl
@@ -5,19 +5,23 @@ This library contains electrical components to build up analog circuits.
module Electrical
using ModelingToolkit, Symbolics, IfElse
+using ModelingToolkit: t, D
+using DynamicQuantities
using ..Thermal: HeatPort
using ..Mechanical.Rotational: Flange, Support
using ..Blocks: RealInput, RealOutput
+using ..DynamicQuantities
+using ..DynamicQuantities: @u_str
-@parameters t
-D = Differential(t)
+import ..rad
+import ..S
export Pin, OnePort
include("utils.jl")
export Capacitor,
- Ground, Inductor, Resistor, Conductor, Short, IdealOpAmp, EMF,
- HeatingResistor
+ Ground, Inductor, Resistor, Conductor, Short, IdealOpAmp, EMF,
+ HeatingResistor
include("Analog/ideal_components.jl")
export CurrentSensor, PotentialSensor, VoltageSensor, PowerSensor, MultiSensor
@@ -38,12 +42,12 @@ export Logic
include("Digital/logic.jl")
export StdLogicVector, StdULogicVector,
- std_ulogic, UX01, UX01Z, X01, X01Z,
- get_logic_level
+ std_ulogic, UX01, UX01Z, X01, X01Z,
+ get_logic_level
include("Digital/logic_vectors.jl")
export LogicTable,
- AndTable, OrTable, NotTable, XorTable
+ AndTable, OrTable, NotTable, XorTable
include("Digital/tables.jl")
end
diff --git a/src/Electrical/utils.jl b/src/Electrical/utils.jl
index ba7d0f541..08cdc3f23 100644
--- a/src/Electrical/utils.jl
+++ b/src/Electrical/utils.jl
@@ -1,6 +1,6 @@
@connector Pin begin
- v(t) # Potential at the pin [V]
- i(t), [connect = Flow] # Current flowing into the pin [A]
+ v(t), [description = "Voltage", unit = u"V"] # Potential at the pin [V]
+ i(t), [description = "Current", unit = u"A", connect = Flow] # Current flowing into the pin [A]
end
@doc """
Pin(; name)
@@ -33,8 +33,8 @@ Component with two electrical pins `p` and `n` and current `i` flows from `p` to
n = Pin()
end
@variables begin
- v(t) = 0.0
- i(t) = 0.0
+ v(t) = 0.0, [description = "Voltage", unit = u"V"]
+ i(t) = 0.0, [description = "Current", unit = u"A", connect = Flow]
end
@equations begin
v ~ p.v - n.v
@@ -70,10 +70,10 @@ Current `i1` flows from `p1` to `n1` and `i2` from `p2` to `n2`.
n2 = Pin()
end
@variables begin
- v1(t) = 0.0
- i1(t) = 0.0
- v2(t) = 0.0
- i2(t) = 0.0
+ v1(t) = 0.0, [description = "Voltage", unit = u"V"]
+ i1(t) = 0.0, [description = "Current", unit = u"A"]
+ v2(t) = 0.0, [description = "Voltage", unit = u"V"]
+ i2(t) = 0.0, [description = "Current", unit = u"A"]
end
@equations begin
v1 ~ p1.v - n1.v
@@ -89,7 +89,7 @@ end
@variables val(t) v(t) i(t)
eqs = [
val ~ IfElse.ifelse((0.0 <= v) & (v <= 0.8) | (2.0 <= v) & (v <= 5.0),
- IfElse.ifelse(v > 2.0, 1, 0), X),
+ IfElse.ifelse(v > 2.0, 1, 0), X)
]
ODESystem(Equation[], t, [val, v, i], [], defaults = Dict(val => 0, i => 0),
name = name)
diff --git a/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl b/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl
index 41236e1de..b58647d84 100644
--- a/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl
+++ b/src/Hydraulic/IsothermalCompressible/IsothermalCompressible.jl
@@ -4,15 +4,13 @@ Library to model iso-thermal compressible liquid fluid flow
module IsothermalCompressible
using ModelingToolkit, Symbolics
+using ModelingToolkit: t_nounits as t, D_nounits as D
using ...Blocks: RealInput, RealOutput
using ...Mechanical.Translational: MechanicalPort, Mass
using IfElse: ifelse
-@parameters t
-D = Differential(t)
-
export HydraulicPort, HydraulicFluid
include("utils.jl")
diff --git a/src/Hydraulic/IsothermalCompressible/components.jl b/src/Hydraulic/IsothermalCompressible/components.jl
index d5ded6358..8e41b39fb 100644
--- a/src/Hydraulic/IsothermalCompressible/components.jl
+++ b/src/Hydraulic/IsothermalCompressible/components.jl
@@ -20,7 +20,7 @@ Caps a hydraulic port to prevent mass flow in or out.
end
eqs = [port.p ~ p
- port.dm ~ 0]
+ port.dm ~ 0]
ODESystem(eqs, t, vars, pars; name, systems)
end
@@ -38,7 +38,7 @@ end
end
eqs = [port.p ~ p
- port.dm ~ dm]
+ port.dm ~ dm]
ODESystem(eqs, t, vars, pars; name, systems)
end
@@ -65,9 +65,9 @@ Variable length internal flow model of the fully developed incompressible flow f
- `port_b`: hydraulic port
"""
@component function TubeBase(add_inertia = true, variable_length = true; p_int, area,
- length_int, head_factor = 1,
- perimeter = 2 * sqrt(area * pi),
- shape_factor = 64, name)
+ length_int, head_factor = 1,
+ perimeter = 2 * sqrt(area * pi),
+ shape_factor = 64, name)
pars = @parameters begin
p_int = p_int
area = area
@@ -116,7 +116,8 @@ Variable length internal flow model of the fully developed incompressible flow f
0
end
- eqs = [0 ~ port_a.dm + port_b.dm]
+ eqs = [0 ~ port_a.dm + port_b.dm
+ domain_connect(port_a, port_b)]
if variable_length
push!(eqs, Δp ~ ifelse(c > 0, shear + inertia, zero(c)))
@@ -149,8 +150,8 @@ Constant length internal flow model discretized by `N` (`FixedVolume`: `N`, `Tub
- `port_b`: hydraulic port
"""
@component function Tube(N, add_inertia = true; p_int, area, length, head_factor = 1,
- perimeter = 2 * sqrt(area * pi),
- shape_factor = 64, name)
+ perimeter = 2 * sqrt(area * pi),
+ shape_factor = 64, name)
@assert(N>0,
"the Tube component must be defined with at least 1 segment (i.e. N>0), found N=$N")
@@ -203,7 +204,7 @@ Constant length internal flow model discretized by `N` (`FixedVolume`: `N`, `Tub
end
eqs = [connect(volumes[1].port, pipe_bases[1].port_a, port_a)
- connect(volumes[end].port, pipe_bases[end].port_b, port_b)]
+ connect(volumes[end].port, pipe_bases[end].port_b, port_b)]
for i in 2:(N - 1)
push!(eqs,
@@ -252,17 +253,17 @@ Reduces the flow from `port_a` to `port_b` by `n`. Useful for modeling parallel
end
eqs = [connect(port_a, port_b, open.port)
- dm_a ~ port_a.dm
- dm_b ~ dm_a / n
- open.dm ~ dm_a - dm_b # extra flow dumps into an open port
- # port_b.dm ~ dm_b # divided flow goes to port_b
- ]
+ dm_a ~ port_a.dm
+ dm_b ~ dm_a / n
+ open.dm ~ dm_a - dm_b # extra flow dumps into an open port
+ # port_b.dm ~ dm_b # divided flow goes to port_b
+ ]
ODESystem(eqs, t, vars, pars; name, systems)
end
@component function ValveBase(reversible = false; p_a_int, p_b_int, minimum_area = 0,
- area_int, Cd, Cd_reverse = Cd, name)
+ area_int, Cd, Cd_reverse = Cd, name)
pars = @parameters begin
p_a_int = p_a_int
p_b_int = p_b_int
@@ -302,8 +303,9 @@ end
end
eqs = [0 ~ port_a.dm + port_b.dm
- dm ~ regRoot(2 * Δp * ρ / c) * x
- y ~ x]
+ domain_connect(port_a, port_b)
+ dm ~ regRoot(2 * Δp * ρ / c) * x
+ y ~ x]
ODESystem(eqs, t, vars, pars; name, systems)
end
@@ -327,9 +329,9 @@ Valve with `area` input and discharge coefficient `Cd` defined by https://en.wik
- `area`: real input setting the valve `area`. When `reversible = true`, negative input reverses flow direction, otherwise a floor of `minimum_area` is enforced.
"""
@component function Valve(reversible = false; p_a_int, p_b_int,
- area_int, Cd, Cd_reverse = Cd,
- minimum_area = 0,
- name)
+ area_int, Cd, Cd_reverse = Cd,
+ minimum_area = 0,
+ name)
pars = @parameters begin
p_a_int = p_a_int
p_b_int = p_b_int
@@ -350,14 +352,14 @@ Valve with `area` input and discharge coefficient `Cd` defined by https://en.wik
vars = []
eqs = [connect(base.port_a, port_a)
- connect(base.port_b, port_b)
- base.area ~ area.u]
+ connect(base.port_b, port_b)
+ base.area ~ area.u]
ODESystem(eqs, t, vars, pars; name, systems, defaults = [area.u => area_int])
end
@component function VolumeBase(; p_int, x_int = 0, area, dead_volume = 0, Χ1 = 1, Χ2 = 1,
- name)
+ name)
pars = @parameters begin
p_int = p_int
x_int = x_int
@@ -382,10 +384,10 @@ end
p = port.p
eqs = [vol ~ dead_volume + area * x
- D(x) ~ dx
- D(rho) ~ drho
- rho ~ full_density(port, p)
- dm ~ drho * vol * Χ1 + rho * area * dx * Χ2]
+ D(x) ~ dx
+ D(rho) ~ drho
+ rho ~ full_density(port, p)
+ dm ~ drho * vol * Χ1 + rho * area * dx * Χ2]
ODESystem(eqs, t, vars, pars; name, systems)
end
@@ -422,8 +424,8 @@ Fixed fluid volume.
p = port.p
eqs = [D(rho) ~ drho
- rho ~ full_density(port, p)
- dm ~ drho * vol]
+ rho ~ full_density(port, p)
+ dm ~ drho * vol]
ODESystem(eqs, t, vars, pars; name, systems)
end
@@ -475,27 +477,27 @@ dm ────► │ │ area
- `flange`: mechanical translational port
"""
@component function DynamicVolume(N, add_inertia = true, reversible = false;
- p_int,
- area,
- x_int = 0,
- x_max,
- x_min = 0,
- x_damp = x_min,
- direction = +1,
-
- # Tube
- perimeter = 2 * sqrt(area * pi),
- shape_factor = 64,
- head_factor = 1,
-
- # Valve
- Cd = 1e2,
- Cd_reverse = Cd,
- minimum_area = 0,
- name)
+ p_int,
+ area,
+ x_int = 0,
+ x_max,
+ x_min = 0,
+ x_damp = x_min,
+ direction = +1,
+
+ # Tube
+ perimeter = 2 * sqrt(area * pi),
+ shape_factor = 64,
+ head_factor = 1,
+
+ # Valve
+ Cd = 1e2,
+ Cd_reverse = Cd,
+ minimum_area = 0,
+ name)
@assert(N>=0,
"the Tube component must be defined with 0 or more segments (i.e. N>=0), found N=$N")
- @assert (direction == +1)||(direction == -1) "direction arument must be +/-1, found $direction"
+ @assert (direction == +1)||(direction == -1) "direction argument must be +/-1, found $direction"
#TODO: How to set an assert effective_length >= length ??
pars = @parameters begin
@@ -549,7 +551,7 @@ dm ────► │ │ area
p_int,
x_int = 0,
area,
- dead_volume = N == 0 ? area * x_max : 0,
+ dead_volume = N == 0 ? area * x_int : 0,
Χ1 = N == 0 ? 1 : 0,
Χ2 = 1)
@@ -562,9 +564,9 @@ dm ────► │ │ area
end
eqs = [vol ~ x * area
- D(x) ~ flange.v * direction
- damper.area ~ damper_area
- connect(port, damper.port_b)]
+ D(x) ~ flange.v * direction
+ damper.area ~ damper_area
+ connect(port, damper.port_b)]
volumes = []
if N > 0
@@ -595,8 +597,9 @@ dm ────► │ │ area
for i in 1:N
push!(eqs,
- volumes[i].dx ~ ifelse((vol >= (i - 1) * (x_max / N) * area) &
- (vol < (i) * (x_max / N) * area),
+ volumes[i].dx ~ ifelse(
+ (vol >= (i - 1) * (x_max / N) * area) &
+ (vol < (i) * (x_max / N) * area),
flange.v * direction, 0))
push!(eqs, pipe_bases[i].x ~ volumes[i].vol / volumes[i].area)
end
@@ -635,17 +638,17 @@ end
end
eqs = [D(x) ~ dx
- flange.v ~ dx
- flange.f ~ 0 #TODO: model flow force
- connect(valve.port_a, port_a)
- connect(valve.port_b, port_b)
- valve.area ~ x * 2π * d]
+ flange.v ~ dx
+ flange.f ~ 0 #TODO: model flow force
+ connect(valve.port_a, port_a)
+ connect(valve.port_b, port_b)
+ valve.area ~ x * 2π * d]
ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0])
end
@component function SpoolValve2Way(reversible = false; p_s_int, p_a_int, p_b_int, p_r_int,
- m, g, x_int, Cd, d, name)
+ m, g, x_int, Cd, d, name)
pars = @parameters begin
p_s_int = p_s_int
p_a_int = p_a_int
@@ -679,37 +682,37 @@ end
end
eqs = [connect(vSA.port_a, port_s)
- connect(vSA.port_b, port_a)
- connect(vBR.port_a, port_b)
- connect(vBR.port_b, port_r)
- connect(vSA.flange, vBR.flange, mass.flange, flange)]
+ connect(vSA.port_b, port_a)
+ connect(vBR.port_a, port_b)
+ connect(vBR.port_b, port_r)
+ connect(vSA.flange, vBR.flange, mass.flange, flange)]
ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0])
end
@component function Actuator(N, add_inertia = true, reversible = false;
- p_a_int,
- p_b_int,
- area_a,
- area_b,
- perimeter_a = 2 * sqrt(area_a * pi),
- perimeter_b = 2 * sqrt(area_b * pi),
- length_a_int,
- length_b_int,
- shape_factor_a = 64,
- shape_factor_b = 64,
- head_factor_a = 1,
- head_factor_b = 1,
- m,
- g,
- x_int = 0,
- minimum_volume_a = 0,
- minimum_volume_b = 0,
- damping_volume_a = minimum_volume_a,
- damping_volume_b = minimum_volume_b,
- Cd = 1e4,
- Cd_reverse = Cd,
- name)
+ p_a_int,
+ p_b_int,
+ area_a,
+ area_b,
+ perimeter_a = 2 * sqrt(area_a * pi),
+ perimeter_b = 2 * sqrt(area_b * pi),
+ length_a_int,
+ length_b_int,
+ shape_factor_a = 64,
+ shape_factor_b = 64,
+ head_factor_a = 1,
+ head_factor_b = 1,
+ m,
+ g,
+ x_int = 0,
+ minimum_volume_a = 0,
+ minimum_volume_b = 0,
+ damping_volume_a = minimum_volume_a,
+ damping_volume_b = minimum_volume_b,
+ Cd = 1e4,
+ Cd_reverse = Cd,
+ name)
pars = @parameters begin
p_a_int = p_a_int
p_b_int = p_b_int
@@ -775,10 +778,10 @@ end
end
eqs = [connect(vol_a.port, port_a)
- connect(vol_b.port, port_b)
- connect(vol_a.flange, vol_b.flange, mass.flange, flange)
- D(x) ~ dx
- dx ~ vol_a.flange.v]
+ connect(vol_b.port, port_b)
+ connect(vol_a.flange, vol_b.flange, mass.flange, flange)
+ D(x) ~ dx
+ dx ~ vol_a.flange.v]
ODESystem(eqs, t, vars, pars; name, systems, defaults = [flange.v => 0])
end
diff --git a/src/Hydraulic/IsothermalCompressible/sources.jl b/src/Hydraulic/IsothermalCompressible/sources.jl
index 237962489..4f7f53b68 100644
--- a/src/Hydraulic/IsothermalCompressible/sources.jl
+++ b/src/Hydraulic/IsothermalCompressible/sources.jl
@@ -18,7 +18,7 @@ Hydraulic mass flow input source
vars = []
eqs = [
- port.dm ~ -dm.u,
+ port.dm ~ -dm.u
]
ODESystem(eqs, t, vars, pars; name, systems)
@@ -47,7 +47,7 @@ Fixed pressure source
end
eqs = [
- port.p ~ p,
+ port.p ~ p
]
ODESystem(eqs, t, vars, pars; name, systems)
@@ -79,7 +79,7 @@ input pressure source
end
eqs = [
- port.p ~ p.u,
+ port.p ~ p.u
]
ODESystem(eqs, t, vars, pars; name, systems)
diff --git a/src/Hydraulic/IsothermalCompressible/utils.jl b/src/Hydraulic/IsothermalCompressible/utils.jl
index d3fce9840..bd09012a5 100644
--- a/src/Hydraulic/IsothermalCompressible/utils.jl
+++ b/src/Hydraulic/IsothermalCompressible/utils.jl
@@ -51,8 +51,8 @@ Fluid parameter setter for isothermal compressible fluid domain. Defaults given
- `p_gas`: [Pa] reference pressure (set by `gas_pressure` argument)
"""
@connector function HydraulicFluid(; density = 997, bulk_modulus = 2.09e9,
- viscosity = 0.0010016, gas_density = 0.0073955,
- gas_pressure = -1000, n = 1, let_gas = 1, name)
+ viscosity = 0.0010016, gas_density = 0.0073955,
+ gas_pressure = -1000, n = 1, let_gas = 1, name)
pars = @parameters begin
ρ = density
β = bulk_modulus
@@ -68,7 +68,7 @@ Fluid parameter setter for isothermal compressible fluid domain. Defaults given
end
eqs = [
- dm ~ 0,
+ dm ~ 0
]
ODESystem(eqs, t, vars, pars; name, defaults = [dm => 0])
@@ -119,7 +119,7 @@ end
Symbolics.derivative(::typeof(friction_factor), args, ::Val{1}) = 0
Symbolics.derivative(::typeof(friction_factor), args, ::Val{4}) = 0
function ChainRulesCore.frule(_, ::typeof(friction_factor), args...)
- (friction_factor(args...), ChainRulesCore.ZeroTangent)
+ (friction_factor(args...), ChainRulesCore.ZeroTangent())
end
function transition(x1, x2, y1, y2, x)
diff --git a/src/Magnetic/FluxTubes/FluxTubes.jl b/src/Magnetic/FluxTubes/FluxTubes.jl
index c082b5f0a..698c9dae9 100644
--- a/src/Magnetic/FluxTubes/FluxTubes.jl
+++ b/src/Magnetic/FluxTubes/FluxTubes.jl
@@ -1,15 +1,15 @@
module FluxTubes
using ModelingToolkit
+using ModelingToolkit: t, D
using ...Electrical: Pin
-
-@parameters t
-D = Differential(t)
+import ...Wb
+using ...DynamicQuantities: @u_str
export PositiveMagneticPort, NegativeMagneticPort, TwoPort
include("utils.jl")
export Ground, Idle, Short, Crossing, ConstantPermeance, ConstantReluctance, EddyCurrent,
- ElectroMagneticConverter
+ ElectroMagneticConverter
include("basic.jl")
export ConstantMagneticPotentialDifference, ConstantMagneticFlux
diff --git a/src/Magnetic/FluxTubes/basic.jl b/src/Magnetic/FluxTubes/basic.jl
index e19ab198c..5e673b143 100644
--- a/src/Magnetic/FluxTubes/basic.jl
+++ b/src/Magnetic/FluxTubes/basic.jl
@@ -20,7 +20,7 @@ Idle running branch.
@mtkmodel Idle begin
@extend (Phi,) = two_port = TwoPort()
@equations begin
- Phi ~ 0
+ Phi ~ 0.0
end
end
@@ -68,7 +68,7 @@ Constant permeance.
@mtkmodel ConstantPermeance begin
@extend V_m, Phi = two_port = TwoPort()
@parameters begin
- G_m = 1.0, [description = "Magnetic permeance"]
+ G_m = 1.0, [description = "Magnetic permeance", unit = u"H"]
end
@equations begin
Phi ~ G_m * V_m
@@ -87,7 +87,7 @@ Constant reluctance.
@mtkmodel ConstantReluctance begin
@extend V_m, Phi = two_port = TwoPort(; Phi = 0.0)
@parameters begin
- R_m = 1.0, [description = "Magnetic reluctance"]
+ R_m = 1.0, [description = "Magnetic reluctance", unit = u"H^-1"]
end
@equations begin
V_m ~ Phi * R_m
@@ -114,11 +114,11 @@ Initial magnetic flux flowing into the port_p can be set with `Phi` ([Wb])
N, [description = "Number of turns"]
end
@variables begin
- v(t)
- i(t)
- Phi
+ v(t),
+ [description = "Voltage difference from positive to negative pin", unit = u"V"]
+ i(t), [description = "Current", unit = u"A"]
end
- @extend V_m, Phi = two_port = TwoPort(; Phi = Phi)
+ @extend V_m, Phi = two_port = TwoPort(; Phi)
@components begin
p = Pin()
n = Pin()
@@ -147,15 +147,16 @@ Initial magnetic flux flowing into the port_p can be set with `Phi` ([`Wb`])
"""
@mtkmodel EddyCurrent begin
@variables begin
- Phi
+ Phi, [description = "Magnetic flux", unit = Wb]
end
@parameters begin
- rho = 0.098e-6, [description = "Resistivity of flux tube material"]
- l = 1, [description = "Average length of eddy current path"]
- A = 1, [description = "Cross sectional area of eddy current path"]
- R = rho * l / A # Electrical resistance of eddy current path
+ rho = 0.098e-6, [description = "Resistivity of flux tube material", unit = u"Ω*m"]
+ l = 1, [description = "Average length of eddy current path", unit = u"m"]
+ A = 1, [description = "Cross sectional area of eddy current path", unit = u"m^2"]
+ R = rho * l / A,
+ [description = "Electrical resistance of eddy current path", unit = u"Ω"]
end
- @extend (V_m, Phi) = two_port = TwoPort(; Phi = Phi)
+ @extend (V_m, Phi) = two_port = TwoPort(; Phi)
@equations begin
D(Phi) ~ V_m * R
end
diff --git a/src/Magnetic/FluxTubes/sources.jl b/src/Magnetic/FluxTubes/sources.jl
index d2c5f8dfe..dff4d228d 100644
--- a/src/Magnetic/FluxTubes/sources.jl
+++ b/src/Magnetic/FluxTubes/sources.jl
@@ -7,16 +7,17 @@ Parameters:
- `V_m`: [A] Magnetic potential difference
"""
+
@mtkmodel ConstantMagneticPotentialDifference begin
@components begin
port_p = PositiveMagneticPort()
port_n = NegativeMagneticPort()
end
@parameters begin
- V_m = 0.0, [description = "Magnetic potential difference"]
+ V_m = 0.0, [description = "Magnetic potential difference", unit = u"A"]
end
@variables begin
- Phi(t)
+ Phi(t), [description = "Magnetic flux", unit = u"Wb"]
end
@equations begin
V_m ~ port_p.V_m - port_n.V_m
@@ -40,10 +41,10 @@ Parameters:
port_n = NegativeMagneticPort()
end
@parameters begin
- Phi = 0.0, [description = "Magnetic flux"]
+ Phi = 0.0, [description = "Magnetic flux", unit = u"Wb"]
end
@variables begin
- V_m(t)
+ V_m(t), [description = "Magnetic potential difference", unit = u"A"]
end
@equations begin
V_m ~ port_p.V_m - port_n.V_m
diff --git a/src/Magnetic/FluxTubes/utils.jl b/src/Magnetic/FluxTubes/utils.jl
index f6154d81f..15ee132d6 100644
--- a/src/Magnetic/FluxTubes/utils.jl
+++ b/src/Magnetic/FluxTubes/utils.jl
@@ -1,6 +1,7 @@
@connector MagneticPort begin
- V_m(t), [description = "Magnetic potential at the port"]
- Phi(t), [connect = Flow, description = "Magnetic flux flowing into the port"]
+ V_m(t), [description = "Magnetic potential at the port", unit = u"A"]
+ Phi(t),
+ [connect = Flow, description = "Magnetic flux flowing into the port", unit = u"Wb"]
end
Base.@doc "Port for a Magnetic system." MagneticPort
@@ -30,8 +31,12 @@ Partial component with magnetic potential difference between two magnetic ports
port_n = NegativeMagneticPort()
end
@variables begin
- V_m(t) = 0.0
- Phi(t) = 0.0
+ V_m(t) = 0.0, [description = "Magnetic potential at the port", unit = u"A"]
+ Phi(t) = 0.0, [
+ connect = Flow,
+ description = "Magnetic flux flowing into the port",
+ unit = u"Wb",
+ ]
end
@equations begin
V_m ~ port_p.V_m - port_n.V_m
diff --git a/src/Mechanical/MultiBody2D/MultiBody2D.jl b/src/Mechanical/MultiBody2D/MultiBody2D.jl
index ed9ccb80d..9a0c048b8 100644
--- a/src/Mechanical/MultiBody2D/MultiBody2D.jl
+++ b/src/Mechanical/MultiBody2D/MultiBody2D.jl
@@ -1,10 +1,10 @@
module MultiBody2D
using ModelingToolkit, Symbolics, IfElse
-using ..Translational
+using ModelingToolkit: t, D
+using DynamicQuantities
-@parameters t
-D = Differential(t)
+using ..TranslationalPosition
export Link
include("components.jl")
diff --git a/src/Mechanical/MultiBody2D/components.jl b/src/Mechanical/MultiBody2D/components.jl
index 99fadd123..3a67dfd5c 100644
--- a/src/Mechanical/MultiBody2D/components.jl
+++ b/src/Mechanical/MultiBody2D/components.jl
@@ -41,11 +41,11 @@
end
@components begin
- TX1 = MechanicalPort()
- TY1 = MechanicalPort()
+ TX1 = Flange()
+ TY1 = Flange()
- TX2 = MechanicalPort()
- TY2 = MechanicalPort()
+ TX2 = Flange()
+ TY2 = Flange()
end
@equations begin
@@ -76,12 +76,12 @@
x_cm ~ l * cos(A) / 2 + x1
y_cm ~ l * sin(A) / 2 + y1
TX1.f ~ fx1
- TX1.v ~ dx1
+ TX1.s ~ x1
TY1.f ~ fy1
- TY1.v ~ dy1
+ TY1.s ~ y1
TX2.f ~ fx2
- TX2.v ~ dx2
+ TX2.s ~ x2
TY2.f ~ fy2
- TY2.v ~ dy2
+ TY2.s ~ y2
end
end
diff --git a/src/Mechanical/Rotational/Rotational.jl b/src/Mechanical/Rotational/Rotational.jl
index 9a6fda31f..0c172a9fb 100644
--- a/src/Mechanical/Rotational/Rotational.jl
+++ b/src/Mechanical/Rotational/Rotational.jl
@@ -4,11 +4,14 @@ Library to model 1-dimensional, rotational mechanical systems
module Rotational
using ModelingToolkit, Symbolics, IfElse
+using ModelingToolkit: t, D
+using DynamicQuantities
+
using ...Blocks: RealInput, RealOutput
import ...@symcheck
-
-@parameters t
-D = Differential(t)
+using ...DynamicQuantities: @u_str
+import ...Wb
+import ...rad
export Flange, Support
include("utils.jl")
@@ -22,4 +25,5 @@ include("sources.jl")
export AngleSensor, SpeedSensor, TorqueSensor, RelSpeedSensor
include("sensors.jl")
+
end
diff --git a/src/Mechanical/Rotational/components.jl b/src/Mechanical/Rotational/components.jl
index 6df5d3b1c..ed80aea50 100644
--- a/src/Mechanical/Rotational/components.jl
+++ b/src/Mechanical/Rotational/components.jl
@@ -16,7 +16,7 @@ Flange fixed in housing at a given angle.
flange = Flange()
end
@parameters begin
- phi0 = 0.0, [description = "Fixed offset angle of flange"]
+ phi0 = 0.0, [description = "Fixed offset angle of flange", unit = u"rad"]
end
@equations begin
flange.phi ~ phi0
@@ -45,7 +45,7 @@ end
"""
@mtkmodel Inertia begin
@parameters begin
- J, [description = "Moment of inertia"]
+ J, [description = "Moment of inertia", unit = u"kg*m^2"]
end
@components begin
flange_a = Flange()
@@ -55,9 +55,9 @@ end
@symcheck J > 0 || throw(ArgumentError("Expected `J` to be positive"))
end
@variables begin
- phi(t) = 0.0, [description = "Absolute rotation angle"]
- w(t) = 0.0, [description = "Absolute angular velocity"]
- a(t) = 0.0, [description = "Absolute angular acceleration"]
+ phi(t) = 0.0, [description = "Absolute rotation angle", unit = u"rad"]
+ w(t) = 0.0, [description = "Absolute angular velocity", unit = u"rad*s^-1"]
+ a(t) = 0.0, [description = "Absolute angular acceleration", unit = u"rad*s^-2"]
end
@equations begin
phi ~ flange_a.phi
@@ -94,8 +94,8 @@ Linear 1D rotational spring
@symcheck c > 0 || throw(ArgumentError("Expected `c` to be positive"))
end
@parameters begin
- c, [description = "Spring constant"]
- phi_rel0 = 0.0, [description = "Unstretched spring angle"]
+ c, [description = "Spring constant", unit = u"N*m*rad^-1"]
+ phi_rel0 = 0.0, [description = "Unstretched spring angle", unit = u"rad"]
end
@equations begin
tau ~ c * (phi_rel - phi_rel0)
@@ -129,7 +129,7 @@ Linear 1D rotational damper
@symcheck d > 0 || throw(ArgumentError("Expected `d` to be positive"))
end
@parameters begin
- d, [description = "Damping constant"]
+ d, [description = "Damping constant", unit = u"N*m*s*rad^-1"]
end
@equations begin
tau ~ d * w_rel
@@ -159,10 +159,10 @@ Linear 1D rotational spring and damper
- `phi_rel0`: [`rad`] Unstretched spring angle
"""
@mtkmodel SpringDamper begin
- @extend phi_rel, w_rel, tau = partial_comp = PartialCompliantWithRelativeStates()
+ @extend phi_rel, w_rel, a_rel, tau = partial_comp = PartialCompliantWithRelativeStates()
@variables begin
- tau_c(t), [description = "Spring torque"]
- tau_d(t), [description = "Damper torque"]
+ tau_c(t), [description = "Spring torque", unit = u"N*m"]
+ tau_d(t), [description = "Damper torque", unit = u"N*m"]
end
@parameters begin
d, [description = "Damping constant"]
@@ -197,19 +197,17 @@ This element characterizes any type of gear box which is fixed in the ground and
# Parameters:
- `ratio`: Transmission ratio (flange_a.phi/flange_b.phi)
- - `use_support`: If support flange enabled, otherwise implicitly grounded
+ - `use_support`: If support flange enabled, otherwise implicitly grounded. By default it is `false`
"""
-@mtkmodel IdealGear begin#(; name, ratio, use_support = false)
- @parameters begin
- use_support
- end
- @extend phi_support, flange_a, flange_b = partial_element = PartialElementaryTwoFlangesAndSupport2(use_support = use_support)
+@mtkmodel IdealGear begin
+ @extend phi_support, flange_a, flange_b = partial_element = PartialElementaryTwoFlangesAndSupport2(;
+ use_support = false)
@parameters begin
ratio, [description = "Transmission ratio"]
end
@variables begin
- phi_a(t) = 0.0, [description = "Relative angle between shaft a and the support"]
- phi_b(t) = 0.0, [description = "Relative angle between shaft b and the support"]
+ phi_a(t) = 0.0, [description = "Relative angle between shaft a and the support", unit = u"rad"]
+ phi_b(t) = 0.0, [description = "Relative angle between shaft b and the support", unit = u"rad"]
end
@equations begin
phi_a ~ flange_a.phi - phi_support
@@ -247,14 +245,13 @@ Friction model: "Armstrong, B. and C.C. de Wit, Friction Modeling and Compensati
- `tau_brk`: [`N⋅m`] Breakaway friction torque
"""
@mtkmodel RotationalFriction begin
- @extend w_rel, tau = partial_comp = PartialCompliantWithRelativeStates()
+ @extend phi_rel, w_rel, a_rel, tau = partial_comp = PartialCompliantWithRelativeStates()
@parameters begin
- f, [description = "Viscous friction coefficient"]
- tau_c, [description = "Coulomb friction torque"]
- w_brk, [description = "Breakaway friction velocity"]
- tau_brk, [description = "Breakaway friction torque"]
+ f, [description = "Viscous friction coefficient", unit = u"N*m*s/rad"]
+ tau_c, [description = "Coulomb friction torque", unit = u"N*m"]
+ w_brk, [description = "Breakaway friction velocity", unit = u"rad*s^-1"]
+ tau_brk, [description = "Breakaway friction torque", unit = "N*m"]
end
-
begin
str_scale = sqrt(2 * exp(1)) * (tau_brk - tau_c)
w_st = w_brk * sqrt(2)
diff --git a/src/Mechanical/Rotational/sensors.jl b/src/Mechanical/Rotational/sensors.jl
index d3be89feb..b90c5271c 100644
--- a/src/Mechanical/Rotational/sensors.jl
+++ b/src/Mechanical/Rotational/sensors.jl
@@ -11,7 +11,7 @@ Ideal sensor to measure the absolute flange angle
@mtkmodel AngleSensor begin
@components begin
flange = Flange()
- phi = RealOutput()
+ phi = RealOutput(unit = u"rad")
end
@equations begin
phi.u ~ flange.phi
@@ -32,7 +32,7 @@ Ideal sensor to measure the absolute flange angular velocity
@mtkmodel SpeedSensor begin
@components begin
flange = Flange()
- w = RealOutput()
+ w = RealOutput(unit = u"rad/s")
end
@equations begin
D(flange.phi) ~ w.u
@@ -55,7 +55,7 @@ Ideal sensor to measure the torque between two flanges (`= flange_a.tau`)
@components begin
flange_a = Flange()
flange_b = Flange()
- tau = RealOutput()
+ tau = RealOutput(unit = u"N*m")
end
@equations begin
flange_a.phi ~ flange_b.phi
@@ -72,16 +72,16 @@ Ideal sensor to measure the relative angular velocity
- `flange_a`: [Flange](@ref) Flange of shaft from which sensor information shall be measured
- `flange_b`: [Flange](@ref) Flange of shaft from which sensor information shall be measured
- - `w`: [RealOutput](@ref) Absolute angular velocity of flange
+ - `w_rel`: [RealOutput](@ref) Relative angular velocity
"""
@mtkmodel RelSpeedSensor begin
@components begin
flange_a = Flange()
flange_b = Flange()
- w_rel = RealOutput()
+ w_rel = RealOutput(unit = u"rad/s")
end
@variables begin
- phi_rel(t) = 0.0
+ phi_rel(t) = 0.0, [description = "Relative rotation angle", unit = u"rad"]
end
@equations begin
0 ~ flange_a.tau + flange_b.tau
diff --git a/src/Mechanical/Rotational/sources.jl b/src/Mechanical/Rotational/sources.jl
index 35abf2c52..a0fd00c67 100644
--- a/src/Mechanical/Rotational/sources.jl
+++ b/src/Mechanical/Rotational/sources.jl
@@ -1,12 +1,11 @@
@mtkmodel PartialTorque begin
- @parameters begin
- use_support
- end
- @extend flange, phi_support = partial_element = PartialElementaryOneFlangeAndSupport2(use_support = use_support)
+ @extend flange, phi_support = partial_element = PartialElementaryOneFlangeAndSupport2(;
+ use_support = false)
@variables begin
phi(t),
[
description = "Angle of flange with respect to support (= flange.phi - support.phi)",
+ unit = u"rad",
]
end
@equations begin
@@ -15,7 +14,7 @@
end
"""
- Torque(; name, use_support)
+ Torque(; name, use_support = false)
Input signal acting as external torque on a flange
@@ -33,12 +32,10 @@ Input signal acting as external torque on a flange
- `use_support`
"""
@mtkmodel Torque begin
- @parameters begin
- use_support
- end
- @extend (flange,) = partial_element = PartialElementaryOneFlangeAndSupport2(use_support = use_support)
+ @extend (flange,) = partial_element = PartialElementaryOneFlangeAndSupport2(;
+ use_support = false)
@components begin
- tau = RealInput()
+ tau = RealInput(unit = u"N*m")
end
@equations begin
flange.tau ~ -tau.u
@@ -46,7 +43,7 @@ Input signal acting as external torque on a flange
end
"""
- ConstantTorque(; name, tau_constant, use_support)
+ ConstantTorque(; name, tau_constant, use_support = false)
Constant torque source
@@ -61,21 +58,28 @@ Constant torque source
# Arguments:
- `tau_constant`: The constant torque applied by the source
-- `use_support`: Whether or not an internal support flange is added.
+- `use_support`: Whether or not an internal support flange is added. By default, it is `false`
"""
-@mtkmodel ConstantTorque begin #(; name, tau_constant, use_support = false)
+@mtkmodel ConstantTorque begin
@parameters begin
tau_constant,
[
description = "Constant torque (if negative, torque is acting as load in positive direction of rotation)",
+ unit = u"N*m",
]
- use_support
end
- @extend flange, phi = partial_element = PartialTorque(; use_support = use_support)
+ @extend flange, phi = partial_element = PartialTorque(; use_support = false)
@variables begin
- tau(t), [description = "Accelerating torque acting at flange (= -flange.tau)"]
+ tau(t),
+ [
+ description = "Accelerating torque acting at flange (= -flange.tau)",
+ unit = u"N*m",
+ ]
w(t),
- [description = "Angular velocity of flange with respect to support (= der(phi))"]
+ [
+ description = "Angular velocity of flange with respect to support",
+ unit = u"rad*s^-1",
+ ]
end
@equations begin
w ~ D(phi)
@@ -85,7 +89,7 @@ Constant torque source
end
"""
- Speed(; name, use_support, exact = false, f_crit = 50)
+ Speed(; name, use_support = false, exact = false, f_crit = 50)
Forced movement of a flange according to a reference angular velocity signal
@@ -108,18 +112,23 @@ Forced movement of a flange according to a reference angular velocity signal
@named partial_element = PartialElementaryOneFlangeAndSupport2(use_support = use_support)
@unpack flange, phi_support = partial_element
@named w_ref = RealInput()
- @variables phi(t)=0.0 w(t)=0.0 a(t)=0.0
+ @variables phi(t) = 0.0 [description = "Angle of flange with respect to support", unit = u"rad"]
+ @variables w(t) = 0.0 [description = "Angular velocity", unit = u"rad*s^-1"]
+ @variables a(t) = 0.0 [description = "Angular acceleration", unit = u"rad*s^-2"]
eqs = [phi ~ flange.phi - phi_support
- D(phi) ~ w]
+ D(phi) ~ w]
if exact
pars = []
push!(eqs, w ~ w_ref.u)
push!(eqs, a ~ 0)
else
- pars = @parameters tau_filt = tau_filt
+ pars = @parameters tau_filt=tau_filt [
+ description = "Time constant of low-pass filter",
+ unit = u"rad*s^-1",
+ ]
push!(eqs, D(w) ~ a)
push!(eqs, a ~ (w_ref.u - w) * tau_filt)
end
- return extend(ODESystem(eqs, t, [phi, w, a], pars; name = name, systems = [w_ref]),
+ return extend(ODESystem(eqs, t, vars, pars; name = name, systems = [w_ref]),
partial_element)
end
diff --git a/src/Mechanical/Rotational/utils.jl b/src/Mechanical/Rotational/utils.jl
index b54bfd06d..b76b21b84 100644
--- a/src/Mechanical/Rotational/utils.jl
+++ b/src/Mechanical/Rotational/utils.jl
@@ -1,6 +1,6 @@
@connector Flange begin
- phi(t), [description = "Rotation angle of flange"]
- tau(t), [connect = Flow, description = "Cut torque in flange"]
+ phi(t), [description = "Rotation angle of flange", unit = u"rad"]
+ tau(t), [connect = Flow, description = "Cut torque in flange", unit = u"N*m"]
end
Base.@doc """
@@ -14,8 +14,8 @@ Base.@doc """
""" Flange
@connector Support begin
- phi(t), [description = "Rotation angle of flange"]
- tau(t), [connect = Flow, description = "Cut torque in flange"]
+ phi(t), [description = "Rotation angle of flange", unit = u"rad"]
+ tau(t), [connect = Flow, description = "Cut torque in flange", unit = u"N*m"]
end
# Base.@doc """
@@ -71,8 +71,8 @@ Partial model for the compliant connection of two rotational 1-dim. shaft flange
flange_b = Flange()
end
@variables begin
- phi_rel(t) = 0.0, [description = "Relative rotation angle between flanges"]
- tau(t) = 0.0, [description = "Torque between flanges"]
+ phi_rel(t)= 0.0, [description = "Relative rotation angle between flanges", unit = u"rad"]
+ tau(t) = 0.0, [description = "Torque between flanges", unit = u"N*m"]
end
@equations begin
phi_rel ~ flange_b.phi - flange_a.phi
@@ -104,10 +104,13 @@ Partial model for the compliant connection of two rotational 1-dim. shaft flange
flange_b = Flange()
end
@variables begin
- phi_rel(t) = 0.0, [description = "Relative rotation angle between flanges"]
- w_rel(t) = 0.0, [description = "Relative angular velocity between flanges"]
- a_rel(t) = 0.0, [description = "Relative angular acceleration between flanges"]
- tau(t) = 0.0, [description = "Torque between flanges"]
+ phi_rel(t) = 0.0, [description = "Relative rotation angle between flanges", unit = u"rad"]
+ w_rel(t) = 0.0, [description = "Relative angular velocity between flanges", unit = u"rad*s^-1"]
+ a_rel(t) = 0.0, [
+ description = "Relative angular acceleration between flanges",
+ unit = u"rad*s^-2",
+ ]
+ tau(t) = 0.0, [description = "Torque between flanges", unit = u"N*m"]
end
@equations begin
phi_rel ~ flange_b.phi - flange_a.phi
@@ -138,11 +141,14 @@ Partial model for a component with one rotational 1-dim. shaft flange and a supp
@component function PartialElementaryOneFlangeAndSupport2(; name, use_support = false)
@named flange = Flange()
sys = [flange]
- @variables phi_support(t)=0.0 [description = "Absolute angle of support flange"]
+ @variables phi_support(t)=0.0 [
+ description = "Absolute angle of support flange",
+ unit = u"rad",
+ ]
if use_support
@named support = Support()
eqs = [support.phi ~ phi_support
- support.tau ~ -flange.tau]
+ support.tau ~ -flange.tau]
push!(sys, support)
else
eqs = [phi_support ~ 0]
@@ -173,11 +179,14 @@ Partial model for a component with two rotational 1-dim. shaft flanges and a sup
@named flange_a = Flange()
@named flange_b = Flange()
sys = [flange_a, flange_b]
- @variables phi_support(t)=0.0 [description = "Absolute angle of support flange"]
+ @variables phi_support(t)=0.0 [
+ description = "Absolute angle of support flange",
+ unit = u"rad",
+ ]
if use_support
@named support = Support()
eqs = [support.phi ~ phi_support
- support.tau ~ -flange_a.tau - flange_b.tau]
+ support.tau ~ -flange_a.tau - flange_b.tau]
push!(sys, support)
else
eqs = [phi_support ~ 0]
diff --git a/src/Mechanical/Translational/Translational.jl b/src/Mechanical/Translational/Translational.jl
index d20961a43..b8aa7dbab 100644
--- a/src/Mechanical/Translational/Translational.jl
+++ b/src/Mechanical/Translational/Translational.jl
@@ -4,13 +4,12 @@ Library to model 1-dimensional, translational mechanical systems
module Translational
using ModelingToolkit, Symbolics
-using ModelingToolkit: getdefault
+using ModelingToolkit: getdefault, t, D
+using DynamicQuantities
using ModelingToolkitStandardLibrary.Blocks: RealInput, RealOutput
using IfElse: ifelse
-
-@parameters t
-D = Differential(t)
+using ...DynamicQuantities: @u_str
export MechanicalPort
include("utils.jl")
diff --git a/src/Mechanical/Translational/components.jl b/src/Mechanical/Translational/components.jl
index 959c2aa5e..a3087d505 100644
--- a/src/Mechanical/Translational/components.jl
+++ b/src/Mechanical/Translational/components.jl
@@ -12,7 +12,7 @@ Use to close a system that has un-connected `MechanicalPort`'s where the force s
flange = MechanicalPort()
end
@variables begin
- f(t) = 0.0
+ f(t) = 0.0, [description = "Force", unit = u"N"]
end
@equations begin
flange.f ~ f
@@ -38,7 +38,7 @@ Fixes a flange position (velocity = 0)
end
"""
- Mass(; name, v_0 = 0.0, m, s_0 = nothing, g = nothing)
+ Mass(; name, v_0 = 0.0, m, s = nothing, g = nothing)
Sliding mass with inertia
@@ -46,13 +46,13 @@ Sliding mass with inertia
- `m`: [kg] mass of sliding body
- `v_0`: [m/s] Initial value of absolute linear velocity of sliding mass (default 0 m/s)
- - `s_0`: [m] (optional) initial value of absolute position of sliding mass
+ - `s`: [m] (optional) initial value of absolute position of sliding mass
- `g`: [m/s²] (optional) gravity field acting on the mass, positive value acts in the positive direction
# States:
- `v`: [m/s] absolute linear velocity of sliding mass
- - `s`: [m] (optional with parameter s_0) absolute position of sliding mass
+ - `s`: [m] (optional with parameter s) absolute position of sliding mass
# Connectors:
@@ -60,21 +60,20 @@ Sliding mass with inertia
"""
@component function Mass(; name, v = 0.0, m, s = nothing, g = nothing)
pars = @parameters begin
- m = m
+ m = m, [description = "Mass of sliding body", unit = u"kg"]
end
@named flange = MechanicalPort(; v = v)
- vars = @variables begin
- v(t) = v
- f(t) = 0
- end
+ @variables v(t) = v [description = "Absolute linear velocity of sliding mass", unit = u"m/s"]
+ @variables f(t) = 0 [description = "Force", unit = u"N"]
+ vars = [v, f]
eqs = [flange.v ~ v
- flange.f ~ f]
+ flange.f ~ f]
# gravity option
if g !== nothing
- @parameters g = g
+ @parameters g = g [description = "Gravitational force", unit = u"N"]
push!(pars, g)
push!(eqs, D(v) ~ f / m + g)
else
@@ -83,7 +82,7 @@ Sliding mass with inertia
# position option
if s !== nothing
- @variables s(t) = s
+ @variables s(t) = s [description = "Absolute position of sliding mass", unit = u"m"]
push!(vars, s)
push!(eqs, D(s) ~ v)
@@ -103,7 +102,7 @@ Linear 1D translational spring
# Parameters:
- `k`: [N/m] Spring constant
- - `delta_s`: initial spring stretch
+ - `delta_s`: [m] initial spring stretch
- `va`: [m/s] Initial value of absolute linear velocity at flange_a (default 0 m/s)
- `v_b_0`: [m/s] Initial value of absolute linear velocity at flange_b (default 0 m/s)
@@ -117,49 +116,45 @@ Linear 1D translational spring
end # default
@component function Spring(::Val{:relative}; name, k, delta_s = 0.0, flange_a__v = 0.0,
- flange_b__v = 0.0)
+ flange_b__v = 0.0)
pars = @parameters begin
- k = k
- end
- vars = @variables begin
- delta_s(t) = delta_s
- f(t) = 0
+ k = k, [description = "Spring constant", unit = u"N/m"]
end
+ @variables delta_s(t) = delta_s [description = "Initial spring stretch", unit = u"m"]
+ @variables f(t) = 0 [description = "Force", unit = u"N"]
@named flange_a = MechanicalPort(; v = flange_a__v)
@named flange_b = MechanicalPort(; v = flange_b__v)
eqs = [D(delta_s) ~ flange_a.v - flange_b.v
- f ~ k * delta_s
- flange_a.f ~ +f
- flange_b.f ~ -f]
- return compose(ODESystem(eqs, t, vars, pars; name = name),
+ f ~ k * delta_s
+ flange_a.f ~ +f
+ flange_b.f ~ -f]
+ return compose(ODESystem(eqs, t, [delta_s, f], pars; name = name),
flange_a,
flange_b) #flange_a.f => +k*delta_s, flange_b.f => -k*delta_s
end
const ABS = Val(:absolute)
@component function Spring(::Val{:absolute}; name, k, sa = 0, sb = 0, flange_a__v = 0.0,
- flange_b__v = 0.0, l = 0)
+ flange_b__v = 0.0, l = 0)
pars = @parameters begin
- k = k
+ k = k, [description = "Spring constant", unit = u"N/m"]
l = l
end
- vars = @variables begin
- sa(t) = sa
- sb(t) = sb
- f(t) = 0
- end
+ @variables sa(t) = sa [unit = u"m"]
+ @variables sb(t) = sb [unit = u"m"]
+ @variables f(t) = 0 [description = "Force", unit = u"N"]
@named flange_a = MechanicalPort(; v = flange_a__v)
@named flange_b = MechanicalPort(; v = flange_b__v)
eqs = [D(sa) ~ flange_a.v
- D(sb) ~ flange_b.v
- f ~ k * (sa - sb - l) #delta_s
- flange_a.f ~ +f
- flange_b.f ~ -f]
- return compose(ODESystem(eqs, t, vars, pars; name = name),
+ D(sb) ~ flange_b.v
+ f ~ k * (sa - sb - l) #delta_s
+ flange_a.f ~ +f
+ flange_b.f ~ -f]
+ return compose(ODESystem(eqs, t, [sa, sb, f], pars; name = name),
flange_a,
flange_b) #, flange_a.f => k * (flange_a__s - flange_b__s - l)
end
@@ -180,11 +175,11 @@ Linear 1D translational damper
"""
@mtkmodel Damper begin
@parameters begin
- d
+ d, [description = "Damping constant", unit = u"N*s/m"]
end
@variables begin
- v(t)
- f(t) = 0.0
+ v(t), [description = "Velocity", unit = u"m/s"]
+ f(t) = 0.0, [description = "Force", unit = u"N"]
end
@components begin
diff --git a/src/Mechanical/Translational/sensors.jl b/src/Mechanical/Translational/sensors.jl
index ec778ee58..1dab8adf5 100644
--- a/src/Mechanical/Translational/sensors.jl
+++ b/src/Mechanical/Translational/sensors.jl
@@ -11,7 +11,7 @@ Linear 1D force input sensor.
@mtkmodel ForceSensor begin
@components begin
flange = MechanicalPort()
- output = RealOutput()
+ output = RealOutput(unit = u"N")
end
@equations begin
@@ -36,11 +36,11 @@ Linear 1D position input sensor.
@mtkmodel PositionSensor begin
@components begin
flange = MechanicalPort()
- output = RealOutput()
+ output = RealOutput(unit = u"m")
end
@variables begin
- s(t) = 0.0
+ s(t) = 0.0, [description = "Absolute position", unit = u"m"]
end
@equations begin
diff --git a/src/Mechanical/Translational/sources.jl b/src/Mechanical/Translational/sources.jl
index 71d4f62f4..508a7ccd7 100644
--- a/src/Mechanical/Translational/sources.jl
+++ b/src/Mechanical/Translational/sources.jl
@@ -11,7 +11,7 @@ Linear 1D force input source
@mtkmodel Force begin
@components begin
flange = MechanicalPort(; v = 0.0)
- f = RealInput()
+ f = RealInput(unit = u"N")
end
@equations begin
@@ -22,35 +22,33 @@ end
"""
Position(solves_force = true; name)
-Linear 1D position input source. Set `solves_force=false` to force input force to 0 (i.e. only the position is given, the respective force needed is already provided elsewhere in the model).
+Linear 1D position input source. Set `solves_force=false` to force input force to 0 (i.e. only the position is given, the respective force needed is already provided elsewhere in the model).
# Connectors:
- `flange`: 1-dim. translational flange
- `s`: real input
"""
-@component function Position(solves_force = true; name)
- vars = []
+@mtkmodel Position begin
+ @structural_parameters begin
+ solves_force = true
+ end
- systems = @named begin
+ @components begin
flange = MechanicalPort(; v = 0)
- s = RealInput()
+ s = RealInput(unit = u"m")
end
- eqs = [
- D(s.u) ~ flange.v,
- ]
-
- !solves_force && push!(eqs, 0 ~ flange.f)
-
- ODESystem(eqs, t, vars, [];
- name, systems)
+ @equations begin
+ D(s.u) ~ flange.v
+ 0 ~ flange.f
+ end
end
"""
Velocity(solves_force = true; name)
-Linear 1D position input source. Set `solves_force=false` to force input force to 0 (i.e. only the velocity is given, the respective force needed is already provided elsewhere in the model).
+Linear 1D position input source. Set `solves_force=false` to force input force to 0 (i.e. only the velocity is given, the respective force needed is already provided elsewhere in the model).
# Connectors:
@@ -60,11 +58,11 @@ Linear 1D position input source. Set `solves_force=false` to force input force
@component function Velocity(solves_force = true; name)
systems = @named begin
flange = MechanicalPort(; v = 0)
- v = RealInput()
+ v = RealInput(unit = u"m/s")
end
eqs = [
- v.u ~ flange.v,
+ v.u ~ flange.v
]
!solves_force && push!(eqs, 0 ~ flange.f)
@@ -75,25 +73,30 @@ end
"""
Acceleration(solves_force = true; name)
-Linear 1D position input source. Set `solves_force=false` to force input force to 0 (i.e. only the acceleration is given, the respective force needed is already provided elsewhere in the model).
+Linear 1D position input source. Set `solves_force=false` to force input force to 0 (i.e. only the acceleration is given, the respective force needed is already provided elsewhere in the model).
# Connectors:
- `flange`: 1-dim. translational flange
- `a`: real input
"""
-@component function Acceleration(solves_force = true; name)
- systems = @named begin
+
+@mtkmodel Acceleration begin#(solves_force = true; name)
+ @structural_parameters begin
+ solves_force = true
+ end
+ @variables begin
+ v(t) = 0, [unit = u"m/s"]
+ end
+ @components begin
flange = MechanicalPort(; v = 0)
- a = RealInput()
+ a = RealInput(unit = u"m/s^2")
+ end
+ @equations begin
+ v ~ flange.v
+ D(v) ~ a.u
+ if solves_force
+ 0 ~ flange.f
+ end
end
-
- vars = @variables v(t) = 0
-
- eqs = [v ~ flange.v
- D(v) ~ a.u]
-
- !solves_force && push!(eqs, 0 ~ flange.f)
-
- ODESystem(eqs, t, vars, []; name, systems)
end
diff --git a/src/Mechanical/Translational/utils.jl b/src/Mechanical/Translational/utils.jl
index 7184d6002..085aa3687 100644
--- a/src/Mechanical/Translational/utils.jl
+++ b/src/Mechanical/Translational/utils.jl
@@ -1,6 +1,6 @@
@connector MechanicalPort begin
- v(t) = 0.0
- f(t) = 0.0, [connect = Flow]
+ v(t) = 0.0, [description = "Velocity of the node", unit = u"m/s"]
+ f(t) = 0.0, [connect = Flow, description = "Force entering the node", unit = u"N"]
end
Base.@doc """
MechanicalPort(;name)
diff --git a/src/Mechanical/TranslationalModelica/TranslationalModelica.jl b/src/Mechanical/TranslationalModelica/TranslationalModelica.jl
index 2aec491b1..041338d7c 100644
--- a/src/Mechanical/TranslationalModelica/TranslationalModelica.jl
+++ b/src/Mechanical/TranslationalModelica/TranslationalModelica.jl
@@ -4,10 +4,10 @@ Library to model 1-dimensional, translational mechanical components.
module TranslationalModelica
using ModelingToolkit, Symbolics, IfElse
-using ...Blocks: RealInput, RealOutput
+using ModelingToolkit: t, D
+using DynamicQuantities
-@parameters t
-D = Differential(t)
+using ...Blocks: RealInput, RealOutput
export Flange
include("utils.jl")
diff --git a/src/Mechanical/TranslationalModelica/components.jl b/src/Mechanical/TranslationalModelica/components.jl
index ad8c92b82..f1ea11817 100644
--- a/src/Mechanical/TranslationalModelica/components.jl
+++ b/src/Mechanical/TranslationalModelica/components.jl
@@ -13,7 +13,7 @@ Flange fixed in housing at a given position.
"""
@mtkmodel Fixed begin
@parameters begin
- s0
+ s0 = 0.0, [description = "Absolute position of sliding mass", unit = u"m"]
end
@components begin
@@ -44,15 +44,14 @@ Sliding mass with inertia
- `flange: 1-dim. translational flange of mass`
"""
@mtkmodel Mass begin
+ @extend flange_a, flange_b, s = pr = PartialRigid(; L = 0.0, s = 0.0)
@parameters begin
- m = 0.0, [description = "Mass of sliding mass [kg]"]
+ m = 0.0, [description = "Mass of sliding mass", unit = u"kg"]
end
@variables begin
- s
- v(t) = 0.0, [description = "Absolute linear velocity of sliding mass [m/s]"]
- a(t) = 0.0, [description = "Absolute linear acceleration of sliding mass [m/s^2]"]
+ v(t) = 0.0, [description = "Absolute linear velocity of sliding mass", unit = u"m/s"]
+ a(t) = 0.0, [description = "Absolute linear acceleration of sliding mass", unit = u"m/s^2"]
end
- @extend flange_a, flange_b, s = pr = PartialRigid(; L = 0.0, s = s)
@equations begin
v ~ D(s)
a ~ D(v)
@@ -78,8 +77,8 @@ Linear 1D translational spring
@mtkmodel Spring begin
@extend flange_a, flange_b, s_rel, f = pc = PartialCompliant()
@parameters begin
- c = 0.0, [description = "Spring constant [N/m]"]
- s_rel0 = 0.0, [description = "Unstretched spring length [m]"]
+ c = 0.0, [description = "Spring constant", unit = u"N/m"]
+ s_rel0 = 0.0, [description = "Unstretched spring length", unit = u"m"]
end
@equations begin
@@ -104,10 +103,10 @@ Linear 1D translational damper
@mtkmodel Damper begin
@extend flange_a, flange_b, v_rel, f = pc = PartialCompliantWithRelativeStates()
@parameters begin
- d = 0.0, [description = "Damping constant [Ns/m]"]
+ d = 0.0, [description = "Damping constant", unit = u"N*s/m"]
end
@variables begin
- lossPower(t) = 0.0, [description = "Power dissipated by the damper [W]"]
+ lossPower(t) = 0.0, [description = "Power dissipated by the damper", unit = u"W"]
end
@equations begin
f ~ d * v_rel
diff --git a/src/Mechanical/TranslationalModelica/sources.jl b/src/Mechanical/TranslationalModelica/sources.jl
index 0110ed2d2..c744ed63f 100644
--- a/src/Mechanical/TranslationalModelica/sources.jl
+++ b/src/Mechanical/TranslationalModelica/sources.jl
@@ -1,15 +1,13 @@
"""
- Force(;name)
+ Force(; name, use_support = false)
Input signal acting as external force on a flange
"""
@mtkmodel Force begin
- @parameters begin
- use_support
- end
- @extend (flange,) = partial_element = PartialElementaryOneFlangeAndSupport2(use_support = use_support)
+ @extend (flange,) = partial_element = PartialElementaryOneFlangeAndSupport2(;
+ use_support = false)
@components begin
- f = RealInput() # Accelerating force acting at flange (= -flange.tau)
+ f = RealInput(unit = u"N") # Accelerating force acting at flange (= -flange.tau)
end
@equations begin
flange.f ~ -f.u
diff --git a/src/Mechanical/TranslationalModelica/utils.jl b/src/Mechanical/TranslationalModelica/utils.jl
index a12016087..273c7c02b 100644
--- a/src/Mechanical/TranslationalModelica/utils.jl
+++ b/src/Mechanical/TranslationalModelica/utils.jl
@@ -1,6 +1,6 @@
@connector Flange begin
- s(t)
- f(t), [connect = Flow]
+ s(t), [description = "Absolute position of sliding mass", unit = u"m"]
+ f(t), [description = "Cut force into the flange", unit = u"N", connect = Flow]
end
Base.@doc """
Flange(;name)
@@ -13,8 +13,8 @@ Base.@doc """
""" Flange
@connector Support begin
- s(t)
- f(t), [connect = Flow]
+ s(t), [description = "Absolute position of sliding mass", unit = u"m"]
+ f(t), [description = "Cut force into the flange", unit = u"N", connect = Flow]
end
Base.@doc """
Support(;name)
@@ -46,8 +46,8 @@ Partial model for the compliant connection of two translational 1-dim. flanges.
@mtkmodel PartialCompliant begin
@extend (flange_a, flange_b) = pt = PartialTwoFlanges()
@variables begin
- s_rel(t) = 0.0, [description = "Relative distance between flanges"]
- f(t) = 0.0, [description = "Force between flanges"]
+ s_rel(t) = 0.0, [description = "Relative distance between flanges", unit = u"m"]
+ f(t) = 0.0, [description = "Force between flanges", unit = u"N"]
end
@equations begin
@@ -71,9 +71,9 @@ Partial model for the compliant connection of two translational 1-dim. flanges.
@mtkmodel PartialCompliantWithRelativeStates begin
@extend flange_a, flange_b = pt = PartialTwoFlanges()
@variables begin
- s_rel(t) = 0.0, [description = "Relative distance between flanges"]
- v_rel(t) = 0.0, [description = "Relative linear velocity))"]
- f(t) = 0.0, [description = "Forces between flanges"]
+ s_rel(t) = 0.0, [description = "Relative distance between flanges", unit = u"m"]
+ v_rel(t) = 0.0, [description = "Relative linear velocity", unit = u"m/s"]
+ f(t) = 0.0, [description = "Forces between flanges", unit = u"N"]
end
@equations begin
@@ -99,9 +99,10 @@ Partial model for a component with one translational 1-dim. shaft flange and a s
"""
function PartialElementaryOneFlangeAndSupport2(; name, use_support = false)
@named flange = Flange()
- @variables s_support(t) [description = "Absolute position of support flange"]
+ @variables s_support(t) [description = "Absolute position of support flange", unit = u"m"]
@variables s(t) [
description = "Distance between flange and support (= flange.s - support.s)",
+ unit = u"m"
]
eqs = [s ~ flange.s - s_support]
if use_support
@@ -130,12 +131,12 @@ Partial model for a component with two translational 1-dim. flanges and a suppor
function PartialElementaryTwoFlangesAndSupport2(; name, use_support = false)
@named flange = Flange()
- @variables s_a(t) [description = "Distance between left flange and support"]
- @variables s_b(t) [description = "Distance between right flange and support"]
- @variables s_support(t) [description = "Absolute position of support flange"]
+ @variables s_a(t) [description = "Distance between left flange and support", unit = u"m"]
+ @variables s_b(t) [description = "Distance between right flange and support", unit = u"m"]
+ @variables s_support(t) [description = "Absolute position of support flange", unit = u"m"]
eqs = [s_a ~ flange_a.s - s_support
- s_b ~ flange_b.s - s_support]
+ s_b ~ flange_b.s - s_support]
if use_support
@named support = Support()
push!(eqs, support.f ~ -flange_a.f - flange_b.f)
@@ -149,10 +150,10 @@ end
@mtkmodel PartialRigid begin
@extend flange_a, flange_b = ptf = PartialTwoFlanges()
@variables begin
- s(t) = 0.0, [description = "Absolute position of center of component"]
+ s(t) = 0.0, [description = "Absolute position of center of component", unit = u"m"]
end
@parameters begin
- L = 0.0, [description = "Length of component, from left flange to right flange"]
+ L = 0.0, [description = "Length of component, from left flange to right flange", unit = u"m"]
end
@equations begin
flange_a.s ~ s - L / 2
diff --git a/src/Mechanical/TranslationalPosition/TranslationalPosition.jl b/src/Mechanical/TranslationalPosition/TranslationalPosition.jl
index 72c91835c..bc20d29eb 100644
--- a/src/Mechanical/TranslationalPosition/TranslationalPosition.jl
+++ b/src/Mechanical/TranslationalPosition/TranslationalPosition.jl
@@ -4,10 +4,10 @@ Library to model 1-dimensional, translational mechanical components.
module TranslationalPosition
using ModelingToolkit, Symbolics, IfElse
-using ...Blocks: RealInput, RealOutput
+using ModelingToolkit: t, D
+using DynamicQuantities
-@parameters t
-D = Differential(t)
+using ...Blocks: RealInput, RealOutput
export Flange
include("utils.jl")
diff --git a/src/Mechanical/TranslationalPosition/components.jl b/src/Mechanical/TranslationalPosition/components.jl
index e1901cb24..5bf0d8876 100644
--- a/src/Mechanical/TranslationalPosition/components.jl
+++ b/src/Mechanical/TranslationalPosition/components.jl
@@ -13,7 +13,7 @@ Flange fixed in housing at a given position.
"""
@mtkmodel Fixed begin
@parameters begin
- s_0
+ s_0 = 0, [description = "Fixed offset position of housing", unit = u"m"]
end
@components begin
flange = Flange(; s = s_0)
@@ -31,13 +31,12 @@ Sliding mass with inertia
# Parameters:
- `m`: [kg] Mass of sliding mass
- - `s_0`: [m] Initial value of absolute position of sliding mass
- - `v_0`: [m/s] Initial value of absolute linear velocity of sliding mass
# States:
- - `s`: [m] Absolute position of sliding mass
- - `v`: [m/s] Absolute linear velocity of sliding mass (= der(s))
+ - `s`: [m] Absolute position of sliding mass. It accepts an initial value, which defaults to 0.0.
+ - `v`: [m/s] Absolute linear velocity of sliding mass (= der(s)). It accepts an initial value, which defaults to 0.0.
+ - `f`: [N] Force. It accepts an initial value, which defaults to 0.0.
# Connectors:
@@ -45,12 +44,12 @@ Sliding mass with inertia
"""
@mtkmodel Mass begin
@parameters begin
- m
+ m, [description = "Mass of sliding mass", unit = u"kg"]
end
@variables begin
- s(t) = 0.0
- v(t) = 0.0
- f(t) = 0.0
+ s(t) = 0.0, [description = "Absolute position of sliding mass", unit = u"m"]
+ v(t) = 0.0, [description = "Absolute linear velocity of sliding mass", unit = u"m*s^-1"]
+ f(t) = 0.0, [description = "Force", unit = u"N"]
end
@components begin
flange = Flange(; s = s)
@@ -65,7 +64,7 @@ end
const REL = Val(:relative)
@component function Spring(::Val{:relative}; name, k, va = 0.0, vb = 0.0,
- delta_s = 0, flange_a__s = 0, flange_b__s = 0)
+ delta_s = 0, flange_a__s = 0, flange_b__s = 0)
pars = @parameters begin
k = k
end
@@ -80,19 +79,22 @@ const REL = Val(:relative)
@named flange_b = Flange()
eqs = [D(flange_a.s) ~ va
- D(flange_b.s) ~ vb
- D(delta_s) ~ va - vb
- f ~ k * delta_s
- flange_a.f ~ +f
- flange_b.f ~ -f]
-
- return compose(ODESystem(eqs, t, vars, pars; name = name,
+ D(flange_b.s) ~ vb
+ D(delta_s) ~ va - vb
+ f ~ k * delta_s
+ flange_a.f ~ +f
+ flange_b.f ~ -f]
+
+ return compose(
+ ODESystem(eqs, t, vars, pars; name = name,
defaults = [
flange_a.s => flange_a__s,
flange_b.s => flange_b__s,
flange_a.f => +delta_s * k,
- flange_b.f => -delta_s * k,
- ]), flange_a, flange_b)
+ flange_b.f => -delta_s * k
+ ]),
+ flange_a,
+ flange_b)
end
const ABS = Val(:absolute)
@@ -105,7 +107,7 @@ Linear 1D translational spring
# Parameters:
- `k`: [N/m] Spring constant
- - `l`: Unstretched spring length
+ - `l`: [m] Unstretched spring length
- `flange_a__s`: [m] Initial value of absolute position of flange_a
- `flange_b__s`: [m] Initial value of absolute position of flange_b
@@ -119,24 +121,22 @@ function Spring(; name, k, flange_a__s = 0, flange_b__s = 0, l = 0)
end #default function
@component function Spring(::Val{:absolute};
- name, k, flange_a__s = 0,
- flange_b__s = 0, l = 0)
+ name, k, flange_a__s = 0,
+ flange_b__s = 0, l = 0)
pars = @parameters begin
- k = k
- l = l
- end
- vars = @variables begin
- f(t) = k * (flange_a__s - flange_b__s - l)
+ k = k, [description = "Spring constant", unit = u"N/m"]
+ l = l, [description = "Unstretched spring length", unit = u"m"]
end
+ @variables f(t) = k * (flange_a__s - flange_b__s - l) [description = "Force", unit = u"N"]
@named flange_a = Flange(; s = flange_a__s, f = k * (flange_a__s - flange_b__s - l))
@named flange_b = Flange(; s = flange_a__s, f = -k * (flange_a__s - flange_b__s - l))
eqs = [
- # delta_s ~ flange_a.s - flange_b.s
- f ~ k * (flange_a.s - flange_b.s - l) #delta_s
- flange_a.f ~ +f
- flange_b.f ~ -f]
+ # delta_s ~ flange_a.s - flange_b.s
+ f ~ k * (flange_a.s - flange_b.s - l) #delta_s
+ flange_a.f ~ +f
+ flange_b.f ~ -f]
return compose(ODESystem(eqs, t, vars, pars; name = name), flange_a, flange_b)
end
@@ -158,12 +158,12 @@ Linear 1D translational damper
"""
@mtkmodel Damper begin
@parameters begin
- d
+ d, [description = "Damping constant", unit = u"N*s/m"]
end
@variables begin
- va(t) = 0.0
- vb(t) = 0.0
- f(t) = +(va - vb) * d
+ va(t) = 0.0, [description = "Velocity of flage a", unit = u"m/s"]
+ vb(t) = 0.0, [description = "Velocity of flage b", unit = u"m/s"]
+ f(t) = +(va - vb) * d, [description = "Force", unit = u"N"]
end
@components begin
diff --git a/src/Mechanical/TranslationalPosition/sources.jl b/src/Mechanical/TranslationalPosition/sources.jl
index 67a94a152..c744ed63f 100644
--- a/src/Mechanical/TranslationalPosition/sources.jl
+++ b/src/Mechanical/TranslationalPosition/sources.jl
@@ -1,15 +1,13 @@
"""
- Force(; name)
+ Force(; name, use_support = false)
Input signal acting as external force on a flange
"""
@mtkmodel Force begin
- @parameters begin
- use_support
- end
- @extend (flange,) = partial_element = PartialElementaryOneFlangeAndSupport2(use_support = use_support)
+ @extend (flange,) = partial_element = PartialElementaryOneFlangeAndSupport2(;
+ use_support = false)
@components begin
- f = RealInput() # Accelerating force acting at flange (= -flange.tau)
+ f = RealInput(unit = u"N") # Accelerating force acting at flange (= -flange.tau)
end
@equations begin
flange.f ~ -f.u
diff --git a/src/Mechanical/TranslationalPosition/utils.jl b/src/Mechanical/TranslationalPosition/utils.jl
index 4e3dd86d9..e133c649e 100644
--- a/src/Mechanical/TranslationalPosition/utils.jl
+++ b/src/Mechanical/TranslationalPosition/utils.jl
@@ -1,6 +1,6 @@
@connector Flange begin
- s(t)
- f(t), [connect = Flow]
+ s(t), [description = "Absolute position of flange", unit = u"m"]
+ f(t), [connect = Flow, description = " Cut force into the flange", unit = u"N"]
end
Base.@doc """
Flange(;name)
@@ -13,8 +13,8 @@ Base.@doc """
""" Flange
@connector Support begin
- s(t)
- f(t), [connect = Flow]
+ s(t), [description = "Absolute position of the support/housing", unit = u"m"]
+ f(t), [connect = Flow, description = "Cut force into the flange", unit = u"N"]
end
Base.@doc """
Support(;name)
@@ -36,16 +36,16 @@ Partial model for the compliant connection of two translational 1-dim. flanges.
- `s_rel`: [m] Relative distance (= flange_b.s - flange_a.s). It accepts an initial value, which defaults to 0.0.
- `f`: [N] Force between flanges (= flange_b.f). It accepts an initial value, which defaults to 0.0.
"""
-@mtkmodel PartialCompliant begin#(; name, s_rel_start = 0.0, f_start = 0.0)
+@mtkmodel PartialCompliant begin
@components begin
flange_a = Flange()
flange_b = Flange()
end
@variables begin
- v_a(t) = 0.0
- v_b(t) = 0.0
- s_rel(t) = 0.0
- f(t) = 0.0
+ v_a(t) = 0.0, [description = "Velocity", unit = u"m/s"]
+ v_b(t) = 0.0, [description = "Velocity", unit = u"m/s"]
+ s_rel(t) = 0.0, [description = "Relative distance ", unit = u"m"]
+ f(t) = 0.0, [description = "Force between flanges", unit = u"N"]
end
@equations begin
D(flange_a.s) ~ v_a
@@ -61,18 +61,9 @@ end
Partial model for the compliant connection of two translational 1-dim. flanges.
-# Parameters:
-
- - `s_rel_start`: [m] Initial relative distance
- - `v_rel_start`: [m/s] Initial relative linear velocity (= der(s_rel))
- - `a_rel_start`: [m/s²] Initial relative linear acceleration (= der(v_rel))
- - `f_start`: [N] Initial force between flanges
+ # States:
-# States:
-
- - `s_rel`: [m] Relative distance (= flange_b.phi - flange_a.phi)
- - `v_rel`: [m/s] Relative linear velocity (= der(s_rel))
- - `a_rel`: [m/s²] Relative linear acceleration (= der(v_rel))
+ - `delta_s`: [m]
- `f`: [N] Force between flanges (= flange_b.f)
"""
@mtkmodel PartialCompliantWithRelativeStates begin
@@ -81,8 +72,8 @@ Partial model for the compliant connection of two translational 1-dim. flanges.
flange_b = Flange()
end
@variables begin
- delta_s(t) = 0.0
- f(t) = 0.0
+ delta_s(t) = 0.0, [unit = u"m"]
+ f(t) = 0.0, [description = "Force between flanges", unit = u"N"]
end
@equations begin
delta_s ~ flange_a.s - flange_b.s
@@ -107,11 +98,12 @@ Partial model for a component with one translational 1-dim. shaft flange and a s
@component function PartialElementaryOneFlangeAndSupport2(; name, use_support = false)
@named flange = Flange()
sys = [flange]
- @variables s_support(t)
+ @variables s_support(t),
+ [description = "Absolute position of support flange", unit = u"m"]
if use_support
@named support = Support()
eqs = [support.s ~ s_support
- support.f ~ -flange.f]
+ support.f ~ -flange.f]
push!(sys, support)
else
eqs = [s_support ~ 0]
@@ -120,13 +112,13 @@ Partial model for a component with one translational 1-dim. shaft flange and a s
end
"""
- PartialElementaryTwoFlangesAndSupport2(;name, use_support=false)
+ PartialElementaryTwoFlangesAndSupport2(; name, use_support = false)
Partial model for a component with two translational 1-dim. flanges and a support used for textual modeling, i.e., for elementary models
# Parameters:
- - `use_support`: If support flange enabled, otherwise implicitly grounded
+ - `use_support`: If support flange enabled, otherwise implicitly grounded. By default it is `false`
# States:
@@ -136,11 +128,11 @@ Partial model for a component with two translational 1-dim. flanges and a suppor
@named flange_a = Flange()
@named flange_b = Flange()
sys = [flange_a, flange_b]
- @variables s_support(t) = 0.0
+ @variables s_support(t) = 0.0 [description = "Absolute position of support flange", unit = u"m"]
if use_support
@named support = Support()
eqs = [support.s ~ s_support
- support.f ~ -flange_a.f - flange_b.f]
+ support.f ~ -flange_a.f - flange_b.f]
push!(sys, support)
else
eqs = [s_support ~ 0]
diff --git a/src/ModelingToolkitStandardLibrary.jl b/src/ModelingToolkitStandardLibrary.jl
index d764f8d81..9b0b6259a 100644
--- a/src/ModelingToolkitStandardLibrary.jl
+++ b/src/ModelingToolkitStandardLibrary.jl
@@ -1,14 +1,15 @@
module ModelingToolkitStandardLibrary
import Symbolics: unwrap
+using DynamicQuantities
"""
@symcheck J > 0 || throw(ArgumentError("Expected `J` to be positive"))
-
+
Omits the check expression if the argument `J` is symbolic.
"""
macro symcheck(ex)
ex.args[1].head === :call ||
- error("Expected an expresion on the form sym > val || error()")
+ error("Expected an expression on the form sym > val || error()")
sym = ex.args[1].args[2]
quote
_issymbolic(x) = !(unwrap(x) isa Real)
@@ -16,6 +17,7 @@ macro symcheck(ex)
end
end
+include("utils.jl")
include("Blocks/Blocks.jl")
include("Mechanical/Mechanical.jl")
include("Thermal/Thermal.jl")
diff --git a/src/Thermal/HeatTransfer/ideal_components.jl b/src/Thermal/HeatTransfer/ideal_components.jl
index f4e0e19c0..dccb76321 100644
--- a/src/Thermal/HeatTransfer/ideal_components.jl
+++ b/src/Thermal/HeatTransfer/ideal_components.jl
@@ -21,11 +21,11 @@ Lumped thermal element storing heat
port = HeatPort()
end
@parameters begin
- C, [description = "Heat capacity of element"]
+ C, [description = "Heat capacity of element", unit = u"J/K"]
end
@variables begin
- T(t) = 273.15 + 20
- der_T(t) = 0.0
+ T(t) = 273.15 + 20, [description = "Temperature of element", unit = u"K"]
+ der_T(t) = 0.0, [description = "Time derivative of temperature", unit = u"K/s"]
end
@equations begin
@@ -56,7 +56,7 @@ see [`Element1D`](@ref)
@mtkmodel ThermalConductor begin
@extend Q_flow, dT = element1d = Element1D()
@parameters begin
- G
+ G, [description = "Constant thermal conductance of material", unit = u"W/K"]
end
@equations begin
Q_flow ~ G * dT
@@ -85,7 +85,7 @@ Lumped thermal element transporting heat without storing it.
@mtkmodel ThermalResistor begin
@extend Q_flow, dT = element1d = Element1D()
@parameters begin
- R
+ R, [description = "Constant thermal resistance of material", unit = u"K/W"]
end
@equations begin
dT ~ R * Q_flow
@@ -114,7 +114,7 @@ Lumped thermal element for heat convection.
@mtkmodel ConvectiveConductor begin
@extend Q_flow, dT = convective_element1d = ConvectiveElement1D()
@parameters begin
- G
+ G, [description = "Convective thermal conductance", unit = u"W/K"]
end
@equations begin
Q_flow ~ G * dT
@@ -143,7 +143,7 @@ Lumped thermal element for heat convection.
@mtkmodel ConvectiveResistor begin
@extend Q_flow, dT = convective_element1d = ConvectiveElement1D()
@parameters begin
- R
+ R, [description = "Constant thermal resistance of material", unit = u"K/W"]
end
@equations begin
dT ~ R * Q_flow
@@ -176,7 +176,7 @@ Lumped thermal element for radiation heat transfer.
@extend Q_flow, dT, port_a, port_b = element1d = Element1D()
@parameters begin
- G
+ G, [description = "Net radiation conductance between two surfaces", unit = "m^2"]
end
@equations begin
Q_flow ~ G * sigma * (port_a.T^4 - port_b.T^4)
@@ -205,7 +205,7 @@ This is a model to collect the heat flows from `m` heatports to one single heatp
port_a = [HeatPort(name = Symbol(:port_a, i)) for i in 1:m]
@named port_b = HeatPort()
eqs = [port_b.Q_flow + sum(k -> k.Q_flow, port_a) ~ 0
- port_b.T ~ port_a[1].T]
+ port_b.T ~ port_a[1].T]
for i in 1:(m - 1)
push!(eqs, port_a[i].T ~ port_a[i + 1].T)
end
diff --git a/src/Thermal/HeatTransfer/sensors.jl b/src/Thermal/HeatTransfer/sensors.jl
index 2282af693..3df90b753 100644
--- a/src/Thermal/HeatTransfer/sensors.jl
+++ b/src/Thermal/HeatTransfer/sensors.jl
@@ -20,7 +20,7 @@ lags are associated with this sensor model.
port = HeatPort()
end
@variables begin
- T(t)
+ T(t), [description = "Absolute temperature", unit = u"K"]
end
@equations begin
T ~ port.T
@@ -51,7 +51,7 @@ output signal in kelvin.
port_b = HeatPort()
end
@variables begin
- T(t)
+ T(t), [description = "Relative temperature", unit = u"K"]
end
@equations begin
T ~ port_a.T - port_b.T
@@ -85,7 +85,7 @@ The output signal is positive, if the heat flows from `port_a` to `port_b`.
port_b = HeatPort()
end
@variables begin
- Q_flow(t)
+ Q_flow(t), [connect = Flow, description = "Heat flow rate", unit = u"W"]
end
@equations begin
port_a.T ~ port_b.T
diff --git a/src/Thermal/HeatTransfer/sources.jl b/src/Thermal/HeatTransfer/sources.jl
index 1ee8d36f7..74863af4e 100644
--- a/src/Thermal/HeatTransfer/sources.jl
+++ b/src/Thermal/HeatTransfer/sources.jl
@@ -19,9 +19,9 @@ the component FixedHeatFlow is connected, if parameter `Q_flow` is positive.
"""
@mtkmodel FixedHeatFlow begin
@parameters begin
- Q_flow = 1.0, [description = "Fixed heat flow rate at port"]
- T_ref = 293.15, [description = "Reference temperature"]
- alpha = 0.0, [description = "Temperature coefficient of heat flow rate"]
+ Q_flow = 1.0, [description = "Fixed heat flow rate at port", unit = u"W"]
+ T_ref = 293.15, [description = "Reference temperature", unit = u"K"]
+ alpha = 0.0, [description = "Temperature coefficient of heat flow rate", unit = u"1/K"]
end
@components begin
port = HeatPort()
@@ -52,7 +52,7 @@ This model defines a fixed temperature `T` at its port in kelvin, i.e., it defin
port = HeatPort()
end
@parameters begin
- T, [description = "Fixed temperature boundary condition"]
+ T, [description = "Fixed temperature boundary condition", unit = u"K"]
end
@equations begin
port.T ~ T
@@ -82,12 +82,12 @@ dependent losses (which are given a reference temperature T_ref).
"""
@mtkmodel PrescribedHeatFlow begin
@parameters begin
- T_ref = 293.15, [description = "Reference temperature"]
- alpha = 0.0, [description = "Temperature coefficient of heat flow rate"]
+ T_ref = 293.15, [description = "Reference temperature", unit = u"K"]
+ alpha = 0.0, [description = "Temperature coefficient of heat flow rate", unit = u"1/K"]
end
@components begin
port = HeatPort()
- Q_flow = RealInput()
+ Q_flow = RealInput(; unit = u"W")
end
@equations begin
port.Q_flow ~ -Q_flow.u * (1 + alpha * (port.T - T_ref))
@@ -111,7 +111,7 @@ the temperature at the specified value.
@mtkmodel PrescribedTemperature begin
@components begin
port = HeatPort()
- T = RealInput()
+ T = RealInput(; unit = u"K")
end
@equations begin
port.T ~ T.u
diff --git a/src/Thermal/Thermal.jl b/src/Thermal/Thermal.jl
index 1a2958fd9..7f7048444 100644
--- a/src/Thermal/Thermal.jl
+++ b/src/Thermal/Thermal.jl
@@ -3,17 +3,17 @@ Library of thermal system components to model heat transfer.
"""
module Thermal
using ModelingToolkit, Symbolics, IfElse
-using ...Blocks: RealInput, RealOutput
+using ModelingToolkit: t, D
+using DynamicQuantities
-@parameters t
-D = Differential(t)
+using ...Blocks: RealInput, RealOutput
export HeatPort, Element1D
include("utils.jl")
export BodyRadiation, ConvectiveConductor, ConvectiveResistor, HeatCapacitor,
- ThermalConductor,
- ThermalResistor, ThermalCollector
+ ThermalConductor,
+ ThermalResistor, ThermalCollector
include("HeatTransfer/ideal_components.jl")
export RelativeTemperatureSensor, HeatFlowSensor, TemperatureSensor
diff --git a/src/Thermal/utils.jl b/src/Thermal/utils.jl
index fe50d4c7d..ba38a9731 100644
--- a/src/Thermal/utils.jl
+++ b/src/Thermal/utils.jl
@@ -1,6 +1,6 @@
@connector HeatPort begin
- T(t) = 273.15 + 20.0
- Q_flow(t) = 0.0, [connect = Flow]
+ T(t) = 273.15 + 20.0, [description = "Temperature", unit = u"K"]
+ Q_flow(t) = 0.0, [connect = Flow, description = "Heat flow rate", unit = u"W"]
end
Base.@doc """
HeatPort(; name, T = 273.15 + 20.0, Q_flow = 0.0)
@@ -35,8 +35,8 @@ flow rate through the element from `port_a` to `port_b`, `Q_flow`.
port_b = HeatPort()
end
@variables begin
- dT(t) = 0.0
- Q_flow(t) = 0.0
+ dT(t) = 0.0, [description = "Temperature difference across the component", unit = u"K"]
+ Q_flow(t) = 0.0, [connect = Flow, description = "Heat flow rate", unit = u"W"]
end
@equations begin
dT ~ port_a.T - port_b.T
@@ -69,8 +69,8 @@ flow rate through the element from `solid` to `fluid`, `Q_flow`.
fluid = HeatPort()
end
@variables begin
- dT(t) = 0.0
- Q_flow(t) = 0.0
+ dT(t) = 0.0, [description = "Temperature difference across the component", unit = u"K"]
+ Q_flow(t) = 0.0, [connect = Flow, description = "Heat flow rate", unit = u"W"]
end
@equations begin
dT ~ solid.T - fluid.T
diff --git a/src/utils.jl b/src/utils.jl
new file mode 100644
index 000000000..c98db13f3
--- /dev/null
+++ b/src/utils.jl
@@ -0,0 +1,6 @@
+# Define units
+
+@register_unit Wb u"m^2*kg*s^-2*A^-1"
+@register_unit H u"kg*m^2*s^−2*A^−2"
+@register_unit rad u"1"
+@register_unit S u"1/Ω"
diff --git a/test/Blocks/continuous.jl b/test/Blocks/continuous.jl
index bdd4739c8..cdf1bc241 100644
--- a/test/Blocks/continuous.jl
+++ b/test/Blocks/continuous.jl
@@ -1,12 +1,18 @@
using ModelingToolkit, ModelingToolkitStandardLibrary, OrdinaryDiffEq
+using ModelingToolkit: t
using ModelingToolkitStandardLibrary.Blocks
using OrdinaryDiffEq: ReturnCode.Success
+using DynamicQuantities: @u_str
+using Test
-@parameters t
+<<<<<<< HEAD
+@parameters t [unit = u"s"]
+=======
+>>>>>>> 645278c (refactor: use default t and D from MTKv9)
#=
Testing strategy:
-The general strategy is to test systems using simple intputs where the solution
+The general strategy is to test systems using simple inputs 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.
=#
@@ -27,9 +33,10 @@ end
@named source = Sine(; frequency = 1)
@named int = Integrator(; k = 1)
@named der = Derivative(; k = 1, T = 0.001)
- @named iosys = ODESystem([
+ @named iosys = ODESystem(
+ [
connect(source.output, der.input),
- connect(der.output, int.input),
+ connect(der.output, int.input)
],
t,
systems = [int, source, der])
@@ -95,7 +102,7 @@ end
@named ss = StateSpace(; A, B, C, D, x = zeros(2))
@named c = Constant(; k = 1)
@named model = ODESystem([
- connect(c.output, ss.input),
+ connect(c.output, ss.input)
],
t,
systems = [ss, c])
@@ -115,7 +122,7 @@ end
y0 = [2]
@named ss = StateSpace(; A, B, C, D, x = zeros(2), u0, y0)
@named model = ODESystem([
- connect(c.output, ss.input),
+ connect(c.output, ss.input)
],
t,
systems = [ss, c])
@@ -136,25 +143,25 @@ Second order demo plant
@component function Plant(; name, x = zeros(2))
@named input = RealInput()
@named output = RealOutput()
- D = Differential(t)
sts = @variables x1(t)=x[1] x2(t)=x[2]
eqs = [D(x1) ~ x2
- D(x2) ~ -x1 - 0.5 * x2 + input.u
- output.u ~ 0.9 * 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
re_val = 2
@named ref = Constant(; k = re_val)
- @named pi_controller = PI(int.k = 1, T = 1)
+ @named pi_controller = PI(k = 1, T = 1)
@named plant = Plant()
@named fb = Feedback()
- @named model = ODESystem([
+ @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),
+ connect(pi_controller.ctr_output, plant.input)
],
t,
systems = [pi_controller, plant, ref, fb])
@@ -172,11 +179,12 @@ end
@named pid_controller = PID(k = 3, Ti = 0.5, Td = 1 / 100)
@named plant = Plant()
@named fb = Feedback()
- @named model = ODESystem([
+ @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),
+ connect(pid_controller.ctr_output, plant.input)
],
t,
systems = [pid_controller, plant, ref, fb])
@@ -189,11 +197,12 @@ end
@testset "PI" begin
@named pid_controller = PID(k = 3, Ti = 0.5, Td = false)
- @named model = ODESystem([
+ @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),
+ connect(pid_controller.ctr_output, plant.input)
],
t,
systems = [pid_controller, plant, ref, fb])
@@ -207,11 +216,12 @@ end
@testset "PD" begin
@named pid_controller = PID(k = 10, Ti = false, Td = 1)
- @named model = ODESystem([
+ @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),
+ connect(pid_controller.ctr_output, plant.input)
],
t,
systems = [pid_controller, plant, ref, fb])
@@ -239,12 +249,13 @@ end
# without anti-windup measure
sol = let
- @named model = ODESystem([
+ @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),
+ connect(sat.output, plant.input)
],
t,
systems = [pi_controller, plant, ref, fb, sat])
@@ -255,12 +266,13 @@ end
# with anti-windup measure
sol_lim = let
- @named model = ODESystem([
+ @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),
+ connect(sat.output, plant.input)
],
t,
systems = [pi_controller_lim, plant, ref, fb, sat])
@@ -284,13 +296,15 @@ end
@testset "LimPID" begin
re_val = 1
@named ref = Constant(; k = re_val)
- @named pid_controller = LimPID(k = 3, Ti = 0.5, Td = 1 / 100, u_max = 1.5, u_min = -1.5,
+ @named pid_controller = LimPID(
+ k = 3, Ti = 0.5, Td = 1 / 100, u_max = 1.5, u_min = -1.5,
Ni = 0.1 / 0.5)
@named plant = Plant()
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(ref.output, pid_controller.reference),
connect(plant.output, pid_controller.measurement),
- connect(pid_controller.ctr_output, plant.input),
+ connect(pid_controller.ctr_output, plant.input)
],
t,
systems = [pid_controller, plant, ref])
@@ -307,10 +321,11 @@ end
@testset "PI" begin
@named pid_controller = LimPID(k = 3, Ti = 0.5, Td = false, u_max = 1.5,
u_min = -1.5, Ni = 0.1 / 0.5)
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(ref.output, pid_controller.reference),
connect(plant.output, pid_controller.measurement),
- connect(pid_controller.ctr_output, plant.input),
+ connect(pid_controller.ctr_output, plant.input)
],
t,
systems = [pid_controller, plant, ref])
@@ -327,10 +342,11 @@ end
@testset "PD" begin
@named pid_controller = LimPID(k = 10, Ti = false, Td = 1, u_max = 1.5,
u_min = -1.5)
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(ref.output, pid_controller.reference),
connect(plant.output, pid_controller.measurement),
- connect(pid_controller.ctr_output, plant.input),
+ connect(pid_controller.ctr_output, plant.input)
],
t,
systems = [pid_controller, plant, ref])
@@ -348,10 +364,11 @@ end
@testset "wp" begin
@named pid_controller = LimPID(k = 3, Ti = 0.5, Td = 1 / 100, u_max = 1.5,
u_min = -1.5, Ni = 0.1 / 0.5, wp = 0, wd = 1)
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(ref.output, pid_controller.reference),
connect(plant.output, pid_controller.measurement),
- connect(pid_controller.ctr_output, plant.input),
+ connect(pid_controller.ctr_output, plant.input)
],
t,
systems = [pid_controller, plant, ref])
@@ -369,10 +386,11 @@ end
@testset "wd" begin
@named pid_controller = LimPID(k = 3, Ti = 0.5, Td = 1 / 100, u_max = 1.5,
u_min = -1.5, Ni = 0.1 / 0.5, wp = 1, wd = 0)
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(ref.output, pid_controller.reference),
connect(plant.output, pid_controller.measurement),
- connect(pid_controller.ctr_output, plant.input),
+ connect(pid_controller.ctr_output, plant.input)
],
t,
systems = [pid_controller, plant, ref])
@@ -391,10 +409,11 @@ end
@testset "PI without AWM" begin
@named pid_controller = LimPID(k = 3, Ti = 0.5, Td = false, u_max = 1.5,
u_min = -1.5, Ni = Inf)
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(ref.output, pid_controller.reference),
connect(plant.output, pid_controller.measurement),
- connect(pid_controller.ctr_output, plant.input),
+ connect(pid_controller.ctr_output, plant.input)
],
t,
systems = [pid_controller, plant, ref])
@@ -408,4 +427,71 @@ end
@test sol[plant.output.u][end]≈re_val atol=1e-3 # zero control error after 100s
@test all(-1.5 .<= sol[pid_controller.ctr_output.u] .<= 1.5) # test limit
end
+
+ @testset "TransferFunction" begin
+ pt1_func(t, k, T) = k * (1 - exp(-t / T)) # Known solution to first-order system
+
+ @named c = Constant(; k = 1)
+ @named pt1 = TransferFunction(b = [1.2], a = [3.14, 1])
+ @named iosys = ODESystem(connect(c.output, pt1.input), t, systems = [pt1, c])
+ sys = structural_simplify(iosys)
+ prob = ODEProblem(sys, Pair[pt1.a_end => 1], (0.0, 100.0))
+ sol = solve(prob, Rodas4())
+ @test sol.retcode == Success
+ @test sol[pt1.output.u]≈pt1_func.(sol.t, 1.2, 3.14) atol=1e-3
+
+ # Test logic for a_end by constructing an integrator
+ @named c = Constant(; k = 1)
+ @named pt1 = TransferFunction(b = [1.2], a = [3.14, 0])
+ @named iosys = ODESystem(connect(c.output, pt1.input), t, systems = [pt1, c])
+ sys = structural_simplify(iosys)
+ prob = ODEProblem(sys, Pair[pt1.a_end => 1], (0.0, 100.0))
+ sol = solve(prob, Rodas4())
+ @test sol.retcode == Success
+ @test sol[pt1.output.u] ≈ sol.t .* (1.2 / 3.14)
+ @test sol[pt1.x[1]] ≈ sol.t .* (1 / 3.14) # Test that scaling of state works properly
+
+ # Test higher order
+
+ 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 pt1 = TransferFunction(b = [w^2], a = [1, 2d * w, w^2])
+ @named iosys = ODESystem(connect(c.output, pt1.input), t, systems = [pt1, c])
+ sys = structural_simplify(iosys)
+ prob = ODEProblem(sys, Pair[pt1.a_end => 1], (0.0, 100.0))
+ sol = solve(prob, Rodas4())
+ @test sol.retcode == Success
+ @test sol[pt1.output.u]≈pt2_func.(sol.t, k, w, d) atol=1e-3
+
+ # test zeros (high-pass version of first test)
+ @named c = Constant(; k = 1)
+ @named pt1 = TransferFunction(b = [1, 0], a = [1, 1])
+ @named iosys = ODESystem(connect(c.output, pt1.input), t, systems = [pt1, c])
+ sys = structural_simplify(iosys)
+ prob = ODEProblem(sys, Pair[pt1.a_end => 1], (0.0, 100.0))
+ sol = solve(prob, Rodas4())
+ @test sol.retcode == Success
+ @test sol[pt1.output.u]≈1 .- pt1_func.(sol.t, 1, 1) atol=1e-3
+ @test sol[pt1.x[1]]≈pt1_func.(sol.t, 1, 1) atol=1e-3 # Test that scaling of state works properly
+
+ # Test with no state
+ @named pt1 = TransferFunction(b = [2.7], a = [pi])
+ @named iosys = ODESystem(connect(c.output, pt1.input), t, systems = [pt1, c])
+ sys = structural_simplify(iosys)
+ prob = ODEProblem(sys, Pair[pt1.a_end => 1], (0.0, 100.0))
+ sol = solve(prob, Rodas4())
+ @test sol.retcode == Success
+ @test all(==(2.7 / pi), sol[pt1.output.u])
+ end
end
diff --git a/test/Blocks/math.jl b/test/Blocks/math.jl
index 1ad2317fc..e200e4e36 100644
--- a/test/Blocks/math.jl
+++ b/test/Blocks/math.jl
@@ -1,19 +1,20 @@
using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkit, OrdinaryDiffEq, Test
using ModelingToolkitStandardLibrary.Blocks: _clamp, _dead_zone
-using ModelingToolkit: inputs, unbound_inputs, bound_inputs
+using ModelingToolkit: inputs, unbound_inputs, bound_inputs, t
using OrdinaryDiffEq: ReturnCode.Success
-
-@parameters t
+using DynamicQuantities: @u_str
@testset "Gain" begin
@named c = Constant(; k = 1)
@named gain = Gain(; k = 1)
@named int = Integrator(; k = 1)
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(c.output, gain.input),
- connect(gain.output, int.input),
- ], t, systems = [int, gain, c])
+ 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))
@@ -30,11 +31,12 @@ end
@named gain = Gain(; k = 1)
@named int = Integrator(; k = 1)
@named fb = Feedback(;)
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(c.output, fb.input1),
connect(fb.input2, int.output),
connect(fb.output, gain.input),
- connect(gain.output, int.input),
+ connect(gain.output, int.input)
],
t,
systems = [int, gain, c, fb])
@@ -53,10 +55,11 @@ end
@named c2 = Sine(; frequency = 1)
@named add = Add(;)
@named int = Integrator(; k = 1)
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(c1.output, add.input1),
connect(c2.output, add.input2),
- connect(add.output, int.input),
+ connect(add.output, int.input)
],
t,
systems = [int, add, c1, c2])
@@ -71,10 +74,11 @@ end
k1 = -1
k2 = 2
@named add = Add(; k1 = k1, k2 = k2)
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(c1.output, add.input1),
connect(c2.output, add.input2),
- connect(add.output, int.input),
+ connect(add.output, int.input)
],
t,
systems = [int, add, c1, c2])
@@ -93,11 +97,12 @@ end
@named c3 = Sine(; frequency = 2)
@named add = Add3(;)
@named int = Integrator(; k = 1)
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(c1.output, add.input1),
connect(c2.output, add.input2),
connect(c3.output, add.input3),
- connect(add.output, int.input),
+ connect(add.output, int.input)
],
t,
systems = [int, add, c1, c2, c3])
@@ -113,11 +118,12 @@ end
k2 = 2
k3 = -pi
@named add = Add3(; k1 = k1, k2 = k2, k3 = k3)
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(c1.output, add.input1),
connect(c2.output, add.input2),
connect(c3.output, add.input3),
- connect(add.output, int.input),
+ connect(add.output, int.input)
],
t,
systems = [int, add, c1, c2, c3])
@@ -136,10 +142,11 @@ end
@named c2 = Sine(; frequency = 1)
@named prod = Product(;)
@named int = Integrator(; k = 1)
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(c1.output, prod.input1),
connect(c2.output, prod.input2),
- connect(prod.output, int.input),
+ connect(prod.output, int.input)
],
t,
systems = [int, prod, c1, c2])
@@ -156,10 +163,11 @@ end
@named c2 = Constant(; k = 2)
@named div = Division(;)
@named int = Integrator(; k = 1)
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(c1.output, div.input1),
connect(c2.output, div.input2),
- connect(div.output, int.input),
+ connect(div.output, int.input)
],
t,
systems = [int, div, c1, c2])
@@ -175,9 +183,10 @@ end
@named c = Sine(; frequency = 1)
@named absb = Abs(;)
@named int = Integrator(; k = 1)
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(c.output, absb.input),
- connect(absb.output, int.input),
+ connect(absb.output, int.input)
],
t,
systems = [int, absb, c])
@@ -215,16 +224,18 @@ end
(Sinh, sinh),
(Cosh, cosh),
(Tanh, tanh),
- (Exp, exp),
+ (Exp, exp)
]
@info "testing $block"
@named source = Sine(frequency = 1, amplitude = 0.5)
@named b = block()
@named int = Integrator()
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(source.output, b.input),
- connect(b.output, int.input),
- ], t, systems = [int, b, source])
+ 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())
@@ -239,10 +250,12 @@ end
@named source = Sine(; frequency = 1, offset = 2, amplitude = 0.5)
@named b = block()
@named int = Integrator()
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(source.output, b.input),
- connect(b.output, int.input),
- ], t, systems = [int, b, source])
+ connect(b.output, int.input)
+ ],
+ t, systems = [int, b, source])
sys = structural_simplify(model)
prob = ODEProblem(sys, Pair[int.x => 0.0, b.input.u => 2.0], (0.0, 1.0))
sol = solve(prob, Rodas4())
@@ -257,10 +270,11 @@ end
@named c2 = Sine(; frequency = 1, offset = 1)
@named b = Atan2(;)
@named int = Integrator(; k = 1)
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(c1.output, b.input1),
connect(c2.output, b.input2),
- connect(b.output, int.input),
+ connect(b.output, int.input)
],
t,
systems = [int, b, c1, c2])
diff --git a/test/Blocks/nonlinear.jl b/test/Blocks/nonlinear.jl
index e13a45047..30130ea17 100644
--- a/test/Blocks/nonlinear.jl
+++ b/test/Blocks/nonlinear.jl
@@ -1,18 +1,19 @@
using ModelingToolkit, OrdinaryDiffEq
+using ModelingToolkit: t
using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkitStandardLibrary.Blocks: _clamp, _dead_zone
using OrdinaryDiffEq: ReturnCode.Success
-
-@parameters t
+using DynamicQuantities: @u_str
@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([
+ @named model = ODESystem(
+ [
connect(c.output, int.input),
- connect(int.output, sat.input),
+ connect(int.output, sat.input)
],
t,
systems = [int, c, sat])
@@ -30,9 +31,10 @@ using OrdinaryDiffEq: ReturnCode.Success
@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([
+ @named iosys = ODESystem(
+ [
connect(source.output, lim.input),
- connect(lim.output, int.input),
+ connect(lim.output, int.input)
],
t,
systems = [source, lim, int])
@@ -57,9 +59,10 @@ end
@named c = Constant(; k = 1)
@named int = Integrator(; k = 1)
@named dz = DeadZone(; u_min = -2, u_max = 1)
- @named model = ODESystem([
+ @named model = ODESystem(
+ [
connect(c.output, int.input),
- connect(int.output, dz.input),
+ connect(int.output, dz.input)
],
t,
systems = [int, c, dz])
@@ -76,9 +79,10 @@ end
@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([
+ @named model = ODESystem(
+ [
connect(source.output, dz.input),
- connect(dz.output, int.input),
+ connect(dz.output, int.input)
],
t,
systems = [int, source, dz])
@@ -102,7 +106,7 @@ end
@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),
+ connect(source.output, rl.input)
],
t,
systems = [source, rl])
diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl
index b2df09461..a0f5d775e 100644
--- a/test/Blocks/sources.jl
+++ b/test/Blocks/sources.jl
@@ -1,18 +1,17 @@
using ModelingToolkit, ModelingToolkitStandardLibrary, OrdinaryDiffEq
+using ModelingToolkit: t, D
using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkitStandardLibrary.Blocks: smooth_sin, smooth_cos, smooth_damped_sin,
- smooth_square, smooth_step, smooth_ramp,
- smooth_triangular, triangular, square
+ smooth_square, smooth_step, smooth_ramp,
+ smooth_triangular, triangular, square
using OrdinaryDiffEq: ReturnCode.Success
-
-@parameters t
-D = Differential(t)
+using DynamicQuantities: @u_str
@testset "Constant" begin
@named src = Constant(k = 2)
@named int = Integrator()
@named iosys = ODESystem([
- connect(src.output, int.input),
+ connect(src.output, int.input)
],
t,
systems = [int, src])
@@ -30,10 +29,11 @@ end
vars = @variables y(t)=1 dy(t)=0 ddy(t)=0
@named src = TimeVaryingFunction(f)
@named int = Integrator()
- @named iosys = ODESystem([y ~ src.output.u
- D(y) ~ dy
- D(dy) ~ ddy
- connect(src.output, int.input)],
+ @named iosys = ODESystem(
+ [y ~ src.output.u
+ D(y) ~ dy
+ D(dy) ~ ddy
+ connect(src.output, int.input)],
t,
systems = [int, src])
sys = structural_simplify(iosys)
@@ -63,7 +63,7 @@ end
@named src = Sine(frequency = frequency, amplitude = amplitude, phase = phase,
offset = offset, start_time = start_time)
@named iosys = ODESystem([
- connect(src.output, int.input),
+ connect(src.output, int.input)
],
t,
systems = [int, src])
@@ -82,7 +82,7 @@ end
start_time = start_time,
smooth = true)
@named smooth_iosys = ODESystem([
- connect(smooth_src.output, int.input),
+ connect(smooth_src.output, int.input)
],
t,
systems = [int, smooth_src])
@@ -92,7 +92,8 @@ end
smooth_sol = solve(smooth_prob, Rodas4())
@test sol.retcode == Success
- @test smooth_sol[smooth_src.output.u]≈smooth_sin.(smooth_sol.t, δ, frequency, amplitude,
+ @test smooth_sol[smooth_src.output.u]≈smooth_sin.(
+ smooth_sol.t, δ, frequency, amplitude,
phase, offset, start_time) atol=1e-3
end
@@ -117,7 +118,7 @@ end
start_time = start_time,
smooth = false)
@named iosys = ODESystem([
- connect(src.output, int.input),
+ connect(src.output, int.input)
],
t,
systems = [int, src])
@@ -135,7 +136,7 @@ end
start_time = start_time,
smooth = true)
@named smooth_iosys = ODESystem([
- connect(smooth_src.output, int.input),
+ connect(smooth_src.output, int.input)
],
t,
systems = [int, smooth_src])
@@ -145,7 +146,8 @@ end
smooth_sol = solve(smooth_prob, Rodas4())
@test smooth_sol.retcode == Success
- @test smooth_sol[smooth_src.output.u]≈smooth_cos.(smooth_sol.t, δ, frequency, amplitude,
+ @test smooth_sol[smooth_src.output.u]≈smooth_cos.(
+ smooth_sol.t, δ, frequency, amplitude,
phase, offset, start_time) atol=1e-3
end
@@ -157,7 +159,7 @@ end
@named src = ContinuousClock(offset = offset, start_time = start_time)
@named int = Integrator()
@named iosys = ODESystem([
- connect(src.output, int.input),
+ connect(src.output, int.input)
],
t,
systems = [int, src])
@@ -183,7 +185,7 @@ end
@named src = Ramp(offset = offset, height = height, duration = duration,
start_time = start_time)
@named iosys = ODESystem([
- connect(src.output, int.input),
+ connect(src.output, int.input)
],
t,
systems = [int, src])
@@ -198,7 +200,7 @@ end
@named smooth_src = Ramp(offset = offset, height = height, duration = duration,
start_time = start_time, smooth = true)
@named smooth_iosys = ODESystem([
- connect(smooth_src.output, int.input),
+ connect(smooth_src.output, int.input)
],
t,
systems = [int, smooth_src])
@@ -221,7 +223,7 @@ end
@named src = Step(offset = offset, height = height, start_time = start_time,
smooth = false)
@named iosys = ODESystem([
- connect(src.output, int.input),
+ connect(src.output, int.input)
],
t,
systems = [int, src])
@@ -238,7 +240,7 @@ end
@named src = Step(offset = offset, height = height, start_time = start_time,
duration = duration, smooth = false)
@named iosys = ODESystem([
- connect(src.output, int.input),
+ connect(src.output, int.input)
],
t,
systems = [int, src])
@@ -254,7 +256,7 @@ end
@named smooth_src = Step(offset = offset, height = height, start_time = start_time,
smooth = true)
@named smooth_iosys = ODESystem([
- connect(smooth_src.output, int.input),
+ connect(smooth_src.output, int.input)
],
t,
systems = [int, smooth_src])
@@ -271,7 +273,7 @@ end
@named smooth_src = Step(offset = offset, height = height, start_time = start_time,
smooth = true, duration = duration)
@named smooth_iosys = ODESystem([
- connect(smooth_src.output, int.input),
+ connect(smooth_src.output, int.input)
],
t,
systems = [int, smooth_src])
@@ -297,7 +299,7 @@ end
@named src = Square(frequency = frequency, amplitude = amplitude,
offset = offset, start_time = start_time)
@named iosys = ODESystem([
- connect(src.output, int.input),
+ connect(src.output, int.input)
],
t,
systems = [int, src])
@@ -312,7 +314,7 @@ end
@named smooth_src = Square(frequency = frequency, amplitude = amplitude,
offset = offset, start_time = start_time, smooth = true)
@named smooth_iosys = ODESystem([
- connect(smooth_src.output, int.input),
+ connect(smooth_src.output, int.input)
],
t,
systems = [int, smooth_src])
@@ -337,7 +339,7 @@ end
@named src = Triangular(frequency = frequency, amplitude = amplitude,
offset = offset, start_time = start_time)
@named iosys = ODESystem([
- connect(src.output, int.input),
+ connect(src.output, int.input)
],
t,
systems = [int, src])
@@ -352,7 +354,7 @@ end
@named smooth_src = Triangular(frequency = frequency, amplitude = amplitude,
offset = offset, start_time = start_time, smooth = true)
@named smooth_iosys = ODESystem([
- connect(smooth_src.output, int.input),
+ connect(smooth_src.output, int.input)
],
t,
systems = [int, smooth_src])
@@ -379,7 +381,7 @@ end
phase = phase, offset = offset, start_time = start_time)
@named int = Integrator()
@named iosys = ODESystem([
- connect(src.output, int.input),
+ connect(src.output, int.input)
],
t,
systems = [int, src])
@@ -395,7 +397,7 @@ end
damping = damping, phase = phase, offset = offset,
start_time = start_time, smooth = true)
@named smooth_iosys = ODESystem([
- connect(smooth_src.output, int.input),
+ connect(smooth_src.output, int.input)
],
t,
systems = [int, smooth_src])
@@ -414,30 +416,65 @@ end
dt = 4e-4
t_end = 10.0
- time = 0:dt:t_end
- x = @. time^2 + 1.0
-
- vars = @variables y(t)=1 dy(t)=0 ddy(t)=0
- @named src = SampledData(Float64)
- @named int = Integrator()
- @named iosys = ODESystem([y ~ src.output.u
- D(y) ~ dy
- D(dy) ~ ddy
- connect(src.output, int.input)],
- t,
- systems = [int, src])
- sys = structural_simplify(iosys)
- s = complete(iosys)
- prob = ODEProblem(sys, [], (0.0, t_end), [s.src.buffer => Parameter(x, dt)];
- tofloat = false)
- prob = remake(prob; p = Parameter.(prob.p))
-
- sol = solve(prob, Rodas4(); initializealg = NoInit())
- @test sol.retcode == Success
- @test sol[src.output.u][1] == 1.0 #check correct initial condition
+ time_span = 0:dt:t_end
+ x = @. time_span^2 + 1.0
+
+ @testset "using Parameter type" begin
+ vars = @variables y(t)=1 dy(t)=0 ddy(t)=0
+ @named src = SampledData(Float64)
+ @named int = Integrator()
+ @named iosys = ODESystem(
+ [y ~ src.output.u
+ D(y) ~ dy
+ D(dy) ~ ddy
+ connect(src.output, int.input)],
+ t,
+ systems = [int, src])
+ sys = structural_simplify(iosys)
+ s = complete(iosys)
+ prob = ODEProblem(sys,
+ [],
+ (0.0, t_end),
+ [s.src.buffer => Parameter(x, dt)];
+ tofloat = false)
+ # prob = remake(prob; p = Parameter.(prob.p)) #<-- no longer needed with ModelingToolkit.jl PR #2231
+
+ sol = solve(prob, Rodas4(); initializealg = NoInit())
+ @test sol.retcode == Success
+ @test sol[src.output.u][1] == 1.0 #check correct initial condition
+
+ @test sol(time)[src.output.u]≈x atol=1e-3
+ @test sol[int.output.u][end]≈1 / 3 * 10^3 + 10.0 atol=1e-3 # closed-form solution to integral
+ @test sol[dy][end]≈2 * time[end] atol=1e-3
+ @test sol[ddy][end]≈2 atol=1e-3
+ end
- @test sol(time)[src.output.u]≈x atol=1e-3
- @test sol[int.output.u][end]≈1 / 3 * 10^3 + 10.0 atol=1e-3 # closed-form solution to integral
- @test sol[dy][end]≈2 * time[end] atol=1e-3
- @test sol[ddy][end]≈2 atol=1e-3
+ @testset "using Vector Based" begin
+ vars = @variables y(t)=1 dy(t)=0 ddy(t)=0
+ @named src = SampledData(dt)
+ @named int = Integrator()
+ @named iosys = ODESystem(
+ [y ~ src.output.u
+ D(y) ~ dy
+ D(dy) ~ ddy
+ connect(src.output, int.input)],
+ t,
+ systems = [int, src])
+ sys = structural_simplify(iosys)
+ s = complete(iosys)
+ prob = ODEProblem(sys,
+ [],
+ (0.0, t_end),
+ [s.src.buffer => x, s.src.sample_time => dt];
+ tofloat = false)
+
+ sol = solve(prob, Rodas4(); initializealg = NoInit())
+ @test sol.retcode == Success
+ @test sol[src.output.u][1] == 1.0 #check correct initial condition
+
+ @test sol(time)[src.output.u]≈x atol=1e-3
+ @test sol[int.output.u][end]≈1 / 3 * 10^3 + 10.0 atol=1e-3 # closed-form solution to integral
+ @test sol[dy][end]≈2 * time[end] atol=1e-3
+ @test sol[ddy][end]≈2 atol=1e-3
+ end
end
diff --git a/test/Blocks/test_analysis_points.jl b/test/Blocks/test_analysis_points.jl
index 39cc222ef..79bfd0b62 100644
--- a/test/Blocks/test_analysis_points.jl
+++ b/test/Blocks/test_analysis_points.jl
@@ -2,8 +2,9 @@ using Test, LinearAlgebra
using ModelingToolkit
using ModelingToolkitStandardLibrary.Blocks
using OrdinaryDiffEq
-using ModelingToolkit: get_eqs, vars, @set!, get_iv
+using ModelingToolkit: get_eqs, vars, @set!, get_iv, t
using ControlSystemsBase
+using DynamicQuantities: @u_str
@named P = FirstOrder(k = 1, T = 1)
@named C = Gain(; k = -1)
@@ -14,14 +15,14 @@ t = ModelingToolkit.get_iv(P)
# Test with explicitly created AnalysisPoint
ap = AnalysisPoint(:plant_input)
eqs = [connect(P.output, C.input)
- connect(C.output, ap, P.input)]
+ connect(C.output, ap, P.input)]
sys = ODESystem(eqs, t, systems = [P, C], name = :hej)
ssys = structural_simplify(sys)
prob = ODEProblem(ssys, [P.x => 1], (0, 10))
sol = solve(prob, Rodas5())
-@test norm(sol[1]) >= 1
-@test norm(sol[end]) < 1e-6 # This fails without the feedback through C
+@test norm(sol.u[1]) >= 1
+@test norm(sol.u[end]) < 1e-6 # This fails without the feedback through C
# plot(sol)
matrices, _ = get_sensitivity(sys, ap)
@@ -45,7 +46,7 @@ T = comp_sensitivity(P, C) # or feedback(P*C)
# Test with automatically created analysis point
eqs = [connect(P.output, C.input)
- connect(C.output, :plant_input, P.input)]
+ connect(C.output, :plant_input, P.input)]
sys = ODESystem(eqs, t, systems = [P, C], name = :hej)
matrices, _ = get_sensitivity(sys, :plant_input)
@@ -84,7 +85,7 @@ matrices, _ = linearize(open_sys, [u], [y])
# Test with more than one AnalysisPoint
eqs = [connect(P.output, :plant_output, C.input)
- connect(C.output, :plant_input, P.input)]
+ connect(C.output, :plant_input, P.input)]
sys = ODESystem(eqs, t, systems = [P, C], name = :hej)
matrices, _ = get_sensitivity(sys, :plant_input)
@@ -111,8 +112,8 @@ matrices2, _ = linearize(sys, :plant_input, [P.output.u])
t = ModelingToolkit.get_iv(P)
eqs = [connect(P.output, :plant_output, add.input2)
- connect(add.output, C.input)
- connect(C.output, :plant_input, P.input)]
+ connect(add.output, C.input)
+ connect(C.output, :plant_input, P.input)]
# eqs = [connect(P.output, add.input2)
# connect(add.output, C.input)
@@ -124,7 +125,7 @@ sys_inner = ODESystem(eqs, t, systems = [P, C, add], name = :inner)
@named F = FirstOrder(k = 1, T = 3)
eqs = [connect(r.output, F.input)
- connect(F.output, sys_inner.add.input1)]
+ connect(F.output, sys_inner.add.input1)]
sys_outer = ODESystem(eqs, t, systems = [F, sys_inner, r], name = :outer)
# test first that the structural_simplify works correctly
@@ -143,7 +144,7 @@ lsys = sminreal(ss(matrices...))
matrices_So, _ = get_sensitivity(sys_outer, :inner_plant_output)
lsyso = sminreal(ss(matrices_So...))
-@test lsys == lsyso || lsys == -1 * lsyso * (-1) # Output and input sensitivites are equal for SISO systems
+@test lsys == lsyso || lsys == -1 * lsyso * (-1) # Output and input sensitivities are equal for SISO systems
## A more complicated test case
using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra
@@ -159,12 +160,12 @@ c = 10 # Damping coefficient
@named inertia2 = Inertia(; J = m2)
@named spring = Spring(; c = k)
@named damper = Damper(; d = c)
-@named torque = Torque(; use_support = false)
+@named torque = Torque()
function SystemModel(u = nothing; name = :model)
eqs = [connect(torque.flange, inertia1.flange_a)
- connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
- connect(inertia2.flange_a, spring.flange_b, damper.flange_b)]
+ connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
+ connect(inertia2.flange_a, spring.flange_b, damper.flange_b)]
if u !== nothing
push!(eqs, connect(torque.tau, u.output))
return @named model = ODESystem(eqs, t;
@@ -174,7 +175,7 @@ function SystemModel(u = nothing; name = :model)
inertia2,
spring,
damper,
- u,
+ u
])
end
ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper], name)
@@ -183,7 +184,7 @@ function AngleSensor(; name)
@named flange = Flange()
@named phi = RealOutput()
eqs = [phi.u ~ flange.phi
- flange.tau ~ 0]
+ flange.tau ~ 0]
return ODESystem(eqs, t, [], []; name = name, systems = [flange, phi])
end
@@ -195,17 +196,17 @@ model = SystemModel()
@named er = Add(k2 = -1)
connections = [connect(r.output, :r, filt.input)
- connect(filt.output, er.input1)
- connect(pid.ctr_output, :u, model.torque.tau)
- connect(model.inertia2.flange_b, sensor.flange)
- connect(sensor.phi, :y, er.input2)
- connect(er.output, :e, pid.err_input)]
+ connect(filt.output, er.input1)
+ connect(pid.ctr_output, :u, model.torque.tau)
+ connect(model.inertia2.flange_b, sensor.flange)
+ connect(sensor.phi, :y, er.input2)
+ connect(er.output, :e, pid.err_input)]
closed_loop = ODESystem(connections, t, systems = [model, pid, filt, sensor, r, er],
name = :closed_loop)
prob = ODEProblem(structural_simplify(closed_loop), Pair[], (0.0, 4.0))
-sol = solve(prob, Rodas4())
+sol = solve(prob, Rodas5P(), reltol = 1e-6, abstol = 1e-9)
# plot(
# plot(sol, vars = [filt.y, model.inertia1.phi, model.inertia2.phi]),
# plot(sol, vars = [pid.ctr_output.u], title = "Control signal"),
@@ -217,7 +218,7 @@ lsys = ss(matrices...) |> sminreal
@test lsys.nx == 8
stepres = ControlSystemsBase.step(c2d(lsys, 0.001), 4)
-@test stepres.y[:]≈sol(0:0.001:4, idxs = model.inertia2.phi) rtol=1e-4
+@test Array(stepres.y[:])≈Array(sol(0:0.001:4, idxs = model.inertia2.phi)) rtol=1e-4
# plot(stepres, plotx=true, ploty=true, size=(800, 1200), leftmargin=5Plots.mm)
# plot!(sol, vars = [model.inertia2.phi], sp=1, l=:dash)
@@ -231,25 +232,30 @@ Si = ss(matrices...)
@test tf(So) ≈ tf(Si)
## A simple multi-level system with loop openings
-@parameters t
@named P_inner = FirstOrder(k = 1, T = 1)
@named feedback = Feedback()
@named ref = Step()
-@named sys_inner = ODESystem([connect(P_inner.output, :y, feedback.input2)
- connect(feedback.output, :u, P_inner.input)
- connect(ref.output, :r, feedback.input1)], t,
+@named sys_inner = ODESystem(
+ [connect(P_inner.output, :y, feedback.input2)
+ connect(feedback.output, :u, P_inner.input)
+ connect(ref.output, :r, feedback.input1)],
+ t,
systems = [P_inner, feedback, ref])
Sinner = sminreal(ss(get_sensitivity(sys_inner, :u)[1]...))
-@named sys_inner = ODESystem([connect(P_inner.output, :y, feedback.input2)
- connect(feedback.output, :u, P_inner.input)], t,
+@named sys_inner = ODESystem(
+ [connect(P_inner.output, :y, feedback.input2)
+ connect(feedback.output, :u, P_inner.input)],
+ t,
systems = [P_inner, feedback])
@named P_outer = FirstOrder(k = rand(), T = rand())
-@named sys_outer = ODESystem([connect(sys_inner.P_inner.output, :y2, P_outer.input)
- connect(P_outer.output, :u2, sys_inner.feedback.input1)], t,
+@named sys_outer = ODESystem(
+ [connect(sys_inner.P_inner.output, :y2, P_outer.input)
+ connect(P_outer.output, :u2, sys_inner.feedback.input1)],
+ t,
systems = [P_outer, sys_inner])
Souter = sminreal(ss(get_sensitivity(sys_outer, :sys_inner_u)[1]...))
@@ -278,7 +284,7 @@ D = [0.0 0.0; 0.0 0.0]
Kss = CS.ss(A, B, C, D)
eqs = [connect(P.output, :plant_output, K.input)
- connect(K.output, :plant_input, P.input)]
+ connect(K.output, :plant_input, P.input)]
sys = ODESystem(eqs, t, systems = [P, K], name = :hej)
matrices, _ = Blocks.get_sensitivity(sys, :plant_input)
@@ -308,8 +314,8 @@ G = CS.feedback(Pss, Kss, pos_feedback = true)
t = ModelingToolkit.get_iv(P)
eqs = [connect(P.output, :plant_output, add.input2)
- connect(add.output, C.input)
- connect(C.output, :plant_input, P.input)]
+ connect(add.output, C.input)
+ connect(C.output, :plant_input, P.input)]
sys_inner = ODESystem(eqs, t, systems = [P, C, add], name = :inner)
@@ -317,7 +323,7 @@ sys_inner = ODESystem(eqs, t, systems = [P, C, add], name = :inner)
@named F = FirstOrder(k = 1, T = 3)
eqs = [connect(r.output, F.input)
- connect(F.output, sys_inner.add.input1)]
+ connect(F.output, sys_inner.add.input1)]
sys_outer = ODESystem(eqs, t, systems = [F, sys_inner, r], name = :outer)
matrices, _ = get_sensitivity(sys_outer, [:inner_plant_input, :inner_plant_output])
diff --git a/test/Electrical/analog.jl b/test/Electrical/analog.jl
index 35fce60d5..0182e4d2b 100644
--- a/test/Electrical/analog.jl
+++ b/test/Electrical/analog.jl
@@ -1,16 +1,19 @@
using ModelingToolkitStandardLibrary.Electrical, ModelingToolkit, OrdinaryDiffEq, Test
+using ModelingToolkit: t
using ModelingToolkitStandardLibrary.Blocks: Step,
- Constant, Sine, Cosine, ExpSine, Ramp,
- Square, Triangular
+ Constant, Sine, Cosine, ExpSine, Ramp,
+ Square, Triangular
using ModelingToolkitStandardLibrary.Blocks: square, triangular
using OrdinaryDiffEq: ReturnCode.Success
+using DynamicQuantities: @u_str
-# using Plots
+@parameters t [unit = u"s"]
+D = Differential(t)
-@parameters t
+# using Plots
@testset "sensors" begin
- @named source = Sine(offset = 1, amplitude = 10, frequency = 5)
+ @named source = Sine(offset = 1, amplitude = 10, frequency = 5, output__unit = u"V")
@named voltage = Voltage()
@named resistor = Resistor(R = 1)
@named capacitor = Capacitor(C = 1, v = 0.0)
@@ -21,15 +24,15 @@ using OrdinaryDiffEq: ReturnCode.Success
@named power_sensor = PowerSensor()
connections = [connect(source.output, voltage.V)
- connect(voltage.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, voltage.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)]
+ connect(voltage.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, voltage.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 = [
@@ -40,10 +43,10 @@ using OrdinaryDiffEq: ReturnCode.Success
ground,
voltage_sensor,
current_sensor,
- power_sensor,
+ power_sensor
])
sys = structural_simplify(model)
- prob = ODAEProblem(sys, [], (0.0, 10.0))
+ prob = ODEProblem(sys, [], (0.0, 10.0))
sol = solve(prob, Tsit5())
# Plots.plot(sol; vars=[capacitor.v, voltage_sensor.v])
@@ -57,7 +60,7 @@ end
# simple voltage divider
@testset "voltage divider with a short branch" begin
- @named source = Constant(k = 10)
+ @named source = Constant(k = 10, output.unit = u"V")
@named voltage = Voltage()
@named R0 = Resistor(R = 1e3)
@named R1 = Resistor(R = 1e3)
@@ -66,10 +69,10 @@ end
@named short = Short()
connections = [connect(source.output, voltage.V)
- connect(voltage.p, R1.p)
- connect(R1.n, short.p, R0.p)
- connect(short.n, R0.n, R2.p)
- connect(R2.n, voltage.n, ground.g)]
+ connect(voltage.p, R1.p)
+ connect(R1.n, short.p, R0.p)
+ connect(short.n, R0.n, R2.p)
+ connect(R2.n, voltage.n, ground.g)]
@named model = ODESystem(connections, t,
systems = [R0, R1, R2, source, short, voltage, ground])
@@ -86,21 +89,21 @@ end
# simple RC
@testset "RC" begin
- @named source = Constant(k = 10)
+ @named source = Constant(k = 10, output.unit = u"V")
@named voltage = Voltage()
@named resistor = Resistor(R = 1)
@named capacitor = Capacitor(C = 1, v = 0.0)
@named ground = Ground()
connections = [connect(source.output, voltage.V)
- connect(voltage.p, resistor.p)
- connect(resistor.n, capacitor.p)
- connect(capacitor.n, voltage.n, ground.g)]
+ connect(voltage.p, resistor.p)
+ connect(resistor.n, capacitor.p)
+ connect(capacitor.n, voltage.n, ground.g)]
@named model = ODESystem(connections, t;
systems = [resistor, capacitor, source, voltage, ground])
sys = structural_simplify(model)
- prob = ODAEProblem(sys, Pair[], (0.0, 10.0))
+ prob = ODEProblem(sys, Pair[], (0.0, 10.0))
sol = solve(prob, Tsit5())
# Plots.plot(sol; vars=[source.v, capacitor.v])
@@ -110,21 +113,21 @@ end
# simple RL
@testset "RL" begin
- @named source = Constant(k = 10)
+ @named source = Constant(k = 10, output.unit = u"V")
@named voltage = Voltage()
@named resistor = Resistor(R = 1)
@named inductor = Inductor(L = 1.0, i = 0.0)
@named ground = Ground()
connections = [connect(source.output, voltage.V)
- connect(voltage.p, resistor.p)
- connect(resistor.n, inductor.p)
- connect(inductor.n, voltage.n, ground.g)]
+ connect(voltage.p, resistor.p)
+ connect(resistor.n, inductor.p)
+ connect(inductor.n, voltage.n, ground.g)]
@named model = ODESystem(connections, t;
systems = [resistor, inductor, source, voltage, ground])
sys = structural_simplify(model)
- prob = ODAEProblem(sys, Pair[], (0.0, 10.0))
+ prob = ODEProblem(sys, Pair[], (0.0, 10.0))
sol = solve(prob, Tsit5())
# Plots.plot(sol; vars=[inductor.i, inductor.i])
@@ -135,15 +138,15 @@ end
@testset "RC with voltage sources" begin
R, C = 1, 1
@named voltage = Voltage()
- @named source_const = Constant(k = 10)
+ @named source_const = Constant(k = 10, output.unit = u"V")
@named source_sin = Sine(offset = 1, amplitude = 10, frequency = 2, start_time = 0.5,
- phase = 0)
- @named source_step = Step(offset = 1, height = 10, start_time = 0.5)
+ phase = 0, output__unit = u"V")
+ @named source_step = Step(offset = 1, height = 10, start_time = 0.5, output__unit = u"V")
@named source_tri = Triangular(offset = 1, start_time = 0.5, amplitude = 10,
- frequency = 2)
+ frequency = 2, output__unit = u"V")
@named source_dsin = ExpSine(offset = 1, amplitude = 10, frequency = 2,
- start_time = 0.5, phase = 0, damping = 0.5)
- @named source_ramp = Ramp(offset = 1, height = 10, start_time = 0.5, duration = 1)
+ start_time = 0.5, phase = 0, damping = 0.5, output__unit = u"V")
+ @named source_ramp = Ramp(offset = 1, height = 10, start_time = 0.5, duration = 1, output__unit = u"V")
sources = [source_const, source_sin, source_step, source_tri, source_dsin, source_ramp]
@named resistor = Resistor(; R)
@@ -152,14 +155,14 @@ end
for source in sources
connections = [connect(source.output, voltage.V)
- connect(voltage.p, resistor.p)
- connect(resistor.n, capacitor.p)
- connect(capacitor.n, voltage.n, ground.g)]
+ connect(voltage.p, resistor.p)
+ connect(resistor.n, capacitor.p)
+ connect(capacitor.n, voltage.n, ground.g)]
@named model = ODESystem(connections, t;
systems = [resistor, capacitor, source, ground, voltage])
sys = structural_simplify(model)
- prob = ODAEProblem(sys, Pair[], (0.0, 10.0))
+ prob = ODEProblem(sys, Pair[], (0.0, 10.0))
sol = solve(prob, Tsit5())
@test sol.retcode == Success
sol = solve(prob, Rodas4())
@@ -173,20 +176,20 @@ end
@testset "RC with current sources" begin
start_time = 2
@named current = Current()
- @named source = Step(start_time = 2)
+ @named source = Step(start_time = 2, output__unit = u"A")
@named resistor = Resistor(R = 1)
@named capacitor = Capacitor(C = 1, v = 0.0)
@named ground = Ground()
connections = [connect(source.output, current.I)
- connect(current.p, resistor.n)
- connect(capacitor.n, resistor.p)
- connect(capacitor.p, current.n, ground.g)]
+ connect(current.p, resistor.n)
+ connect(capacitor.n, resistor.p)
+ connect(capacitor.p, current.n, ground.g)]
@named model = ODESystem(connections, t;
systems = [ground, resistor, current, capacitor, source])
sys = structural_simplify(model)
- prob = ODAEProblem(sys, Pair[], (0.0, 10.0))
+ prob = ODEProblem(sys, Pair[], (0.0, 10.0))
sol = solve(prob, Tsit5())
y(x, st) = (x .> st) .* abs.(collect(x) .- st)
@test sol.retcode == Success
@@ -202,17 +205,17 @@ end
@named R2 = Resistor(R = 100 * R)
@named C1 = Capacitor(C = 1 / (2 * pi * f * R), v = 0.0)
@named opamp = IdealOpAmp()
- @named square_source = Square(amplitude = Vin)
+ @named square_source = Square(amplitude = Vin, output__unit = u"V")
@named voltage = Voltage()
@named sensor = VoltageSensor()
connections = [connect(square_source.output, voltage.V)
- connect(voltage.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, voltage.n)
- connect(opamp.p2, sensor.p)
- connect(sensor.n, ground.g)]
+ connect(voltage.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, voltage.n)
+ connect(opamp.p2, sensor.p)
+ connect(sensor.n, ground.g)]
@named model = ODESystem(connections, t,
systems = [
R1,
@@ -222,11 +225,11 @@ end
voltage,
C1,
ground,
- sensor,
+ sensor
])
sys = structural_simplify(model)
u0 = [C1.v => 0.0
- R1.v => 0.0]
+ R1.v => 0.0]
prob = ODEProblem(sys, u0, (0, 100.0))
sol = solve(prob, Rodas4())
@test sol.retcode == Success
@@ -251,26 +254,32 @@ _damped_sine_wave(x, f, A, st, ϕ, d) = exp((st - x) * d) * A * sin(2 * π * f *
@named ground = Ground()
@named voltage = Voltage()
@named voltage_sensor = VoltageSensor()
- @named step = Step(start_time = st, offset = o, height = h)
+ @named step = Step(start_time = st, offset = o, height = h, output__unit = u"V")
@named cosine = Cosine(offset = o, amplitude = A, frequency = f, start_time = st,
- phase = ϕ)
- @named sine = Sine(offset = o, amplitude = A, frequency = f, start_time = st, phase = ϕ)
+ phase = ϕ, output__unit = u"V")
+ @named sine = Sine(offset = o, amplitude = A, frequency = f, start_time = st,
+ phase = ϕ, output__unit = u"V")
@named damped_sine = ExpSine(offset = o, amplitude = A, frequency = f, start_time = st,
- phase = ϕ, damping = d)
- @named ramp = Ramp(offset = o, start_time = st, duration = et - st, height = h)
- @named vsquare = Square(offset = o, start_time = st, amplitude = A, frequency = f)
- @named tri = Triangular(offset = o, start_time = st, amplitude = A, frequency = f)
+ phase = ϕ, damping = d, output__unit = u"V")
+ @named ramp = Ramp(offset = o, start_time = st, duration = et - st, height = h,
+ output__unit = u"V")
+ @named vsquare = Square(offset = o, start_time = st, amplitude = A, frequency = f,
+ output__unit = u"V")
+ @named tri = Triangular(offset = o, start_time = st, amplitude = A, frequency = f,
+ output__unit = u"V")
# @named vsawtooth = SawTooth(amplitude=A, start_time=st, frequency=f, offset=o)
sources = [step, cosine, sine, damped_sine, ramp, tri, vsquare] #, vsawtooth]
function waveforms(i, x)
- getindex([o .+ _step.(x, h, st),
+ getindex(
+ [o .+ _step.(x, h, st),
o .+ (x .> st) .* _cos_wave.(x, f, A, st, ϕ),
o .+ (x .> st) .* _sine_wave.(x, f, A, st, ϕ),
o .+ (x .> st) .* _damped_sine_wave.(x, f, A, st, ϕ, d),
o .+ _ramp.(x, st, (et - st), h),
triangular.(x, f, A, o, st),
- square.(x, f, A, o, st)], i)
+ square.(x, f, A, o, st)],
+ i)
end
# o .+ (x .> st). * _sawtooth_wave.(x, δ, f, A, st),
@@ -278,9 +287,9 @@ _damped_sine_wave(x, f, A, st, ϕ, d) = exp((st - x) * d) * A * sin(2 * π * f *
source = sources[i]
@info "Testing Voltage with $(source.name) source"
eqs = [connect(source.output, voltage.V)
- connect(voltage.p, voltage_sensor.p, res.p)
- connect(res.n, cap.p)
- connect(ground.g, voltage_sensor.n, voltage.n, cap.n)]
+ connect(voltage.p, voltage_sensor.p, res.p)
+ connect(res.n, cap.p)
+ connect(ground.g, voltage_sensor.n, voltage.n, cap.n)]
@named vmodel = ODESystem(eqs, t,
systems = [
voltage_sensor,
@@ -288,13 +297,13 @@ _damped_sine_wave(x, f, A, st, ϕ, d) = exp((st - x) * d) * A * sin(2 * π * f *
cap,
source,
voltage,
- ground,
+ ground
])
vsys = structural_simplify(vmodel)
u0 = [cap.v => 0.0]
- prob = ODAEProblem(vsys, u0, (0, 10.0))
+ prob = ODEProblem(vsys, u0, (0, 10.0))
sol = solve(prob, dt = 0.1, Tsit5())
@test sol.retcode == Success
@@ -314,26 +323,32 @@ end
@named cap = Capacitor(C = 1, v = 0.0)
@named current_sensor = CurrentSensor()
@named current = Current()
- @named step = Step(start_time = st, offset = o, height = h)
+ @named step = Step(start_time = st, offset = o, height = h, output__unit = u"A")
@named cosine = Cosine(offset = o, amplitude = A, frequency = f, start_time = st,
- phase = ϕ)
- @named sine = Sine(offset = o, amplitude = A, frequency = f, start_time = st, phase = ϕ)
+ phase = ϕ, output__unit = u"A")
+ @named sine = Sine(offset = o, amplitude = A, frequency = f, start_time = st,
+ phase = ϕ, output__unit = u"A")
@named damped_sine = ExpSine(offset = o, amplitude = A, frequency = f, start_time = st,
- phase = ϕ, damping = d)
- @named ramp = Ramp(offset = o, start_time = st, duration = et - st, height = h)
- @named vsquare = Square(offset = o, start_time = st, amplitude = A, frequency = f)
- @named tri = Triangular(offset = o, start_time = st, amplitude = A, frequency = f)
+ phase = ϕ, damping = d, output__unit = u"A")
+ @named ramp = Ramp(offset = o, start_time = st, duration = et - st, height = h,
+ output__unit = u"A")
+ @named vsquare = Square(offset = o, start_time = st, amplitude = A, frequency = f,
+ output__unit = u"A")
+ @named tri = Triangular(offset = o, start_time = st, amplitude = A, frequency = f,
+ output__unit = u"A")
# @named isawtooth = SawTooth(amplitude=A, start_time=st, frequency=f, offset=o)
sources = [step, cosine, sine, damped_sine, ramp, tri, vsquare] #, idamped_sine]
function waveforms(i, x)
- getindex([o .+ _step.(x, h, st),
+ getindex(
+ [o .+ _step.(x, h, st),
o .+ (x .> st) .* _cos_wave.(x, f, A, st, ϕ),
o .+ (x .> st) .* _sine_wave.(x, f, A, st, ϕ),
o .+ (x .> st) .* _damped_sine_wave.(x, f, A, st, ϕ, d),
o .+ _ramp.(x, st, (et - st), h),
triangular.(x, f, A, o, st),
- square.(x, f, A, o, st)], i)
+ square.(x, f, A, o, st)],
+ i)
end
# # o .+ (x .> st). * _sawtooth_wave.(x, δ, f, A, st)
@@ -341,10 +356,10 @@ end
source = sources[i]
@info "Testing Current with $(source.name) source"
eqs = [connect(source.output, current.I)
- connect(current.p, current_sensor.n)
- connect(current_sensor.p, res.p)
- connect(res.n, cap.p)
- connect(current.n, ground.g, cap.n)]
+ connect(current.p, current_sensor.n)
+ connect(current_sensor.p, res.p)
+ connect(res.n, cap.p)
+ connect(current.n, ground.g, cap.n)]
@named model = ODESystem(eqs, t,
systems = [
current_sensor,
@@ -352,13 +367,13 @@ end
current,
res,
cap,
- ground,
+ ground
])
isys = structural_simplify(model)
u0 = [cap.v => 0.0]
- prob = ODAEProblem(isys, u0, (0, 10.0))
+ prob = ODEProblem(isys, u0, (0, 10.0))
sol = solve(prob, dt = 0.1, Tsit5())
@test sol.retcode == Success
diff --git a/test/Electrical/digital.jl b/test/Electrical/digital.jl
index 7b5740f4c..0349c5732 100644
--- a/test/Electrical/digital.jl
+++ b/test/Electrical/digital.jl
@@ -4,7 +4,7 @@ using ModelingToolkitStandardLibrary.Electrical: U, X, F0, F1, Z, W, L, H, DC, U
using ModelingToolkitStandardLibrary.Electrical: AndTable, OrTable, NotTable, XorTable
using ModelingToolkitStandardLibrary.Electrical: get_logic_level
using OrdinaryDiffEq: ReturnCode.Success
-
+# using ModelingToolkit: t, D
# using ModelingToolkitStandardLibrary.Electrical: Set, Reset
@testset "Logic, logic-vectors and helpers" begin
@@ -21,7 +21,7 @@ using OrdinaryDiffEq: ReturnCode.Success
@test zero(Logic) == zero(U) == F0
@test one(Logic) == one(U) == F1
@test ones(Logic, 2, 2) == [F1 F1
- F1 F1]
+ F1 F1]
# Logic vectors
u_logic = StdULogicVector([U, W, X, 1])
@@ -32,7 +32,7 @@ using OrdinaryDiffEq: ReturnCode.Success
@test typeof(logic.logic) == Vector{Logic}
@test get_logic_level(logic) == [1, 6, 2, 4]
- # Predefiend logic vectors
+ # Predefined logic vectors
@test std_ulogic.logic == [U, X, F0, F1, Z, W, L, H, DC]
@test UX01.logic == [U, X, F0, F1]
@test UX01Z.logic == [U, X, F0, F1, Z]
@@ -41,7 +41,7 @@ using OrdinaryDiffEq: ReturnCode.Success
# Logic vector helpers
test_logic_matrix = StdULogicVector([U F0
- F1 X])
+ F1 X])
test_logic_vector = StdLogicVector([U, F0, F1, X])
size(test_logic_matrix) == (2, 2)
@@ -96,8 +96,6 @@ end
#=
-@parameters t
-
@named set1 = Set()
@named reset1 = Reset()
@named set2 = Set()
diff --git a/test/Hydraulic/isothermal_compressible.jl b/test/Hydraulic/isothermal_compressible.jl
index 8ccda0303..5e2f7c4ea 100644
--- a/test/Hydraulic/isothermal_compressible.jl
+++ b/test/Hydraulic/isothermal_compressible.jl
@@ -1,13 +1,11 @@
using ModelingToolkit, OrdinaryDiffEq, Test
+using ModelingToolkit: t, D
import ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible as IC
import ModelingToolkitStandardLibrary.Blocks as B
import ModelingToolkitStandardLibrary.Mechanical.Translational as T
using ModelingToolkitStandardLibrary.Blocks: Parameter
-@parameters t
-D = Differential(t)
-
NEWTON = NLNewton(check_div = false, always_new = true, max_iter = 100, relax = 9 // 10)
@testset "Fluid Domain and Tube" begin
@@ -26,9 +24,9 @@ NEWTON = NLNewton(check_div = false, always_new = true, max_iter = 100, relax =
end
eqs = [connect(stp.output, src.p)
- connect(fluid, src.port)
- connect(src.port, res.port_a)
- connect(res.port_b, vol.port)]
+ connect(fluid, src.port)
+ connect(src.port, res.port_a)
+ connect(res.port_b, vol.port)]
ODESystem(eqs, t, [], pars; name, systems)
end
@@ -41,7 +39,7 @@ NEWTON = NLNewton(check_div = false, always_new = true, max_iter = 100, relax =
probs = [ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 0.05))
for sys in syss] #
sols = [solve(prob, ImplicitEuler(nlsolve = NEWTON); initializealg = NoInit(),
- dt = 1e-4, adaptive = false)
+ dt = 1e-4, adaptive = false)
for prob in probs]
s1_2 = complete(sys1_2)
@@ -79,9 +77,9 @@ end
end
eqs = [connect(fluid, sink.port)
- connect(sink.port, valve.port_a)
- connect(valve.port_b, vol.port)
- connect(valve.area, ramp.output)]
+ connect(sink.port, valve.port_a)
+ connect(valve.port_b, vol.port)
+ connect(valve.area, ramp.output)]
ODESystem(eqs, t, [], pars; name, systems)
end
@@ -127,12 +125,12 @@ end
end
eqs = [connect(fluid, src1.port)
- connect(fluid, src2.port)
- connect(src1.port, vol1.port)
- connect(src2.port, vol2.port)
- connect(vol1.flange, mass.flange, vol2.flange)
- connect(src1.p, sin1.output)
- connect(src2.p, sin2.output)]
+ connect(fluid, src2.port)
+ connect(src1.port, vol1.port)
+ connect(src2.port, vol2.port)
+ connect(vol1.flange, mass.flange, vol2.flange)
+ connect(src1.p, sin1.output)
+ connect(src2.p, sin2.output)]
ODESystem(eqs, t, [], pars; name, systems)
end
@@ -195,8 +193,6 @@ end
@testset "Actuator System" begin
function System(use_input, f; name)
- @parameters t
-
pars = @parameters begin
p_s = 200e5
p_r = 5e5
@@ -261,17 +257,17 @@ end
push!(systems, input)
eqs = [connect(input.output, pos.s)
- connect(valve.flange, pos.flange)
- connect(valve.port_a, piston.port_a)
- connect(piston.flange, body.flange)
- connect(piston.port_b, m1.port_a)
- connect(m1.port_b, pipe.port_b)
- connect(pipe.port_a, m2.port_b)
- connect(m2.port_a, valve.port_b)
- connect(src.port, valve.port_s)
- connect(snk.port, valve.port_r)
- connect(fluid, src.port, snk.port)
- D(body.v) ~ ddx]
+ connect(valve.flange, pos.flange)
+ connect(valve.port_a, piston.port_a)
+ connect(piston.flange, body.flange)
+ connect(piston.port_b, m1.port_a)
+ connect(m1.port_b, pipe.port_b)
+ connect(pipe.port_a, m2.port_b)
+ connect(m2.port_a, valve.port_b)
+ connect(src.port, valve.port_s)
+ connect(snk.port, valve.port_r)
+ connect(fluid, src.port, snk.port)
+ D(body.v) ~ ddx]
ODESystem(eqs, t, vars, pars; name, systems)
end
@@ -281,6 +277,12 @@ end
sys = structural_simplify(system)
defs = ModelingToolkit.defaults(sys)
s = complete(system)
+ dt = 1e-4
+ time = 0:dt:0.1
+ x = @. 0.9 * (time > 0.015) * (time - 0.015)^2 - 25 * (time > 0.02) * (time - 0.02)^3
+ defs[s.input.buffer] = Parameter(x, dt)
+ #prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 0.1);
+ # tofloat = false)#, jac = true)
prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 0.1);
tofloat = false, jac = true)
@@ -292,15 +294,9 @@ end
@test Symbol(defs[s.valve.port_r.ρ]) == Symbol(s.fluid.ρ)
@test Symbol(defs[s.snk.port.ρ]) == Symbol(s.fluid.ρ)
- dt = 1e-4
- time = 0:dt:0.1
- x = @. 0.9 * (time > 0.015) * (time - 0.015)^2 - 25 * (time > 0.02) * (time - 0.02)^3
-
- defs[s.input.buffer] = Parameter(x, dt)
# defs[s.piston.Cd_reverse] = 0.1
- p = Parameter.(ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat = false))
- prob = remake(prob; p, tspan = (0, time[end]))
+ prob = remake(prob; tspan = (0, time[end]))
@time sol = solve(prob, ImplicitEuler(nlsolve = NEWTON); adaptive = false, dt,
initializealg = NoInit())
@@ -322,7 +318,7 @@ end
end
eqs = [connect(fluid, cap.port, vol.port)
- connect(vol.flange, mass.flange)]
+ connect(vol.flange, mass.flange)]
ODESystem(eqs, t, [], pars; name, systems)
end
diff --git a/test/Magnetic/magnetic.jl b/test/Magnetic/magnetic.jl
index e378c32f4..673237837 100644
--- a/test/Magnetic/magnetic.jl
+++ b/test/Magnetic/magnetic.jl
@@ -4,11 +4,10 @@ import ModelingToolkitStandardLibrary.Electrical
import ModelingToolkitStandardLibrary.Blocks
import ModelingToolkitStandardLibrary.Magnetic
using ModelingToolkit, OrdinaryDiffEq, Test
+using ModelingToolkit: t, D
using OrdinaryDiffEq: ReturnCode.Success
# using Plots
-@parameters t
-
@testset "Inductor" begin
mu_air = 1
l_air = 0.0001
@@ -26,16 +25,16 @@ using OrdinaryDiffEq: ReturnCode.Success
@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.output, voltage.V)
- connect(voltage.p, r.p)
- connect(r.n, coil.p)
- connect(voltage.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, voltage.n)
- connect(ground_m.port, r_mFe.port_n)]
+ connect(voltage.p, r.p)
+ connect(r.n, coil.p)
+ connect(voltage.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, voltage.n)
+ connect(ground_m.port, r_mFe.port_n)]
@named model = ODESystem(connections, t,
systems = [
source,
@@ -46,7 +45,7 @@ using OrdinaryDiffEq: ReturnCode.Success
r_mAirPar,
r_mFe,
r_mLeak,
- voltage,
+ voltage
])
sys = structural_simplify(model)
prob = ODEProblem(sys, Pair[], (0, 0.1))
diff --git a/test/Mechanical/multibody.jl b/test/Mechanical/multibody.jl
index c00203c84..968c3b03b 100644
--- a/test/Mechanical/multibody.jl
+++ b/test/Mechanical/multibody.jl
@@ -1,28 +1,28 @@
using ModelingToolkit
+using ModelingToolkit: t
using ModelingToolkitStandardLibrary.Mechanical.MultiBody2D
-using ModelingToolkitStandardLibrary.Mechanical.Translational
-using DifferentialEquations
+using ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition
+using OrdinaryDiffEq
# using Setfield
using Test
-@parameters t
-
@named link1 = Link(; m = 1, l = 10, I = 84, g = -9.807)
@named link2 = Link(; m = 1, l = 10, I = 84, g = -9.807, x1_0 = 10)
-@named cart = Mass(; m = 1, s_0 = 0)
+@named cart = Mass(; m = 1, s = 0)
# @named force = SineForce(;amp=3e3, freq=15)
@named fixed = Fixed()
# @named m1 = Mass(;m=0.5)
# @named m2 = Mass(;m=0.5)
eqs = [connect(link1.TX1, cart.flange) #, force.flange)
- connect(link1.TY1, fixed.flange)
- connect(link1.TX2, link2.TX1)
- connect(link1.TY2, link2.TY1)]
+ connect(link1.TY1, fixed.flange)
+ connect(link1.TX2, link2.TX1)
+ connect(link1.TY2, link2.TY1)]
@named model = ODESystem(eqs, t, [], []; systems = [link1, link2, cart, fixed])
sys = structural_simplify(model)
+@test length(unknowns(sys)) == 6
# The below code does work...
#=
diff --git a/test/Mechanical/rotational.jl b/test/Mechanical/rotational.jl
index 79bbbf3ae..5bf33a902 100644
--- a/test/Mechanical/rotational.jl
+++ b/test/Mechanical/rotational.jl
@@ -1,14 +1,12 @@
using ModelingToolkitStandardLibrary.Mechanical.Rotational,
- ModelingToolkit, OrdinaryDiffEq,
- Test
+ ModelingToolkit, OrdinaryDiffEq,
+ Test
+using ModelingToolkit: t, D
import ModelingToolkitStandardLibrary.Blocks
using OrdinaryDiffEq: ReturnCode.Success
# using Plots
-@parameters t
-D = Differential(t)
-
@testset "two inertias" begin
@named fixed = Fixed()
@named inertia1 = Inertia(J = 2) # this one is fixed
@@ -17,8 +15,8 @@ D = Differential(t)
@named inertia2 = Inertia(J = 2, phi = 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)]
+ 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])
@@ -32,7 +30,7 @@ D = Differential(t)
sol = solve(prob, Rodas4())
@test SciMLBase.successful_retcode(sol)
- prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, Pair[], (0, 10.0))
+ prob = DAEProblem(sys, D.(unknowns(sys)) .=> 0.0, Pair[], (0, 10.0))
sol = solve(prob, DFBDF())
@test SciMLBase.successful_retcode(sol)
@test all(sol[inertia1.w] .== 0)
@@ -40,8 +38,8 @@ D = Differential(t)
@named springdamper = SpringDamper(; c = 1e4, d = 10)
connections = [connect(fixed.flange, inertia1.flange_b)
- connect(inertia1.flange_b, springdamper.flange_a)
- connect(springdamper.flange_b, inertia2.flange_a)]
+ connect(inertia1.flange_b, springdamper.flange_a)
+ connect(springdamper.flange_b, inertia2.flange_a)]
@named model = ODESystem(connections, t,
systems = [fixed, inertia1, inertia2, springdamper])
@@ -62,7 +60,7 @@ end
J_motor = 0.1 # Motor inertia
@named fixed = Fixed()
- @named torque = Torque(use_support = true)
+ @named torque = Torque(; use_support = true)
@named inertia1 = Inertia(J = 2, phi = pi / 2)
@named spring = Rotational.Spring(c = 1e4)
@named damper = Damper(d = 10)
@@ -70,10 +68,10 @@ end
@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)]
+ 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 = [
@@ -83,11 +81,11 @@ end
inertia2,
spring,
damper,
- sine,
+ sine
])
sys = structural_simplify(model)
- prob = DAEProblem(sys, D.(states(sys)) .=> 0.0,
- [D(D(inertia2.phi)) => 1.0; D.(states(model)) .=> 0.0], (0, 10.0))
+ prob = DAEProblem(sys, D.(unknowns(sys)) .=> 0.0,
+ [D(D(inertia2.phi)) => 1.0; D.(unknowns(model)) .=> 0.0], (0, 10.0))
sol = solve(prob, DFBDF())
@test SciMLBase.successful_retcode(sol)
@@ -101,9 +99,9 @@ end
## Test with constant torque source
@named torque = ConstantTorque(use_support = true, tau_constant = 1)
connections = [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)]
+ 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 = [
@@ -112,7 +110,7 @@ end
inertia1,
inertia2,
spring,
- damper,
+ damper
])
sys = structural_simplify(model)
@@ -142,15 +140,15 @@ end
@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)]
+ 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 = [
@@ -162,7 +160,7 @@ end
spring,
inertia3,
damper,
- sine,
+ sine
])
@test_skip begin
sys = structural_simplify(model) #key 7 not found
@@ -180,8 +178,8 @@ end
@named lim = Blocks.Limiter(y_max = 6)
@named output = Blocks.RealOutput()
connections = [connect(sine.output, dz.input)
- connect(dz.output, lim.input)
- connect(lim.output, output)]
+ connect(dz.output, lim.input)
+ connect(lim.output, output)]
ODESystem(connections, t, [], []; name = name, systems = [sine, dz, lim, output])
end
@@ -192,15 +190,15 @@ end
@named friction = RotationalFriction(f = 0.001, tau_c = 20, w_brk = 0.06035,
tau_brk = 25)
@named vel_profile = VelocityProfile()
- @named source = Speed(use_support = false)
+ @named source = Speed()
@named angle_sensor = AngleSensor()
connections = [connect(vel_profile.output, source.w_ref)
- connect(source.flange, friction.flange_a)
- connect(friction.flange_b, inertia.flange_a)
- connect(inertia.flange_b, spring.flange_a, damper.flange_a)
- connect(spring.flange_b, damper.flange_b, fixed.flange)
- connect(angle_sensor.flange, inertia.flange_a)]
+ connect(source.flange, friction.flange_a)
+ connect(friction.flange_b, inertia.flange_a)
+ connect(inertia.flange_b, spring.flange_a, damper.flange_a)
+ connect(spring.flange_b, damper.flange_b, fixed.flange)
+ connect(angle_sensor.flange, inertia.flange_a)]
@named model = ODESystem(connections, t,
systems = [
@@ -211,10 +209,10 @@ end
vel_profile,
source,
friction,
- angle_sensor,
+ angle_sensor
])
sys = structural_simplify(model)
- prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, Pair[], (0, 10.0))
+ prob = DAEProblem(sys, D.(unknowns(sys)) .=> 0.0, Pair[], (0, 10.0))
sol = solve(prob, DFBDF())
@test SciMLBase.successful_retcode(sol)
@@ -239,16 +237,16 @@ end
@named rel_speed_sensor = RelSpeedSensor()
connections = [connect(fixed.flange, inertia1.flange_b, rel_speed_sensor.flange_b)
- connect(inertia1.flange_b, torque_sensor.flange_a)
- connect(torque_sensor.flange_b, spring.flange_a, damper.flange_a,
- speed_sensor.flange, rel_speed_sensor.flange_a)
- connect(spring.flange_b, damper.flange_b, inertia2.flange_a)]
+ connect(inertia1.flange_b, torque_sensor.flange_a)
+ connect(torque_sensor.flange_b, spring.flange_a, damper.flange_a,
+ speed_sensor.flange, rel_speed_sensor.flange_a)
+ connect(spring.flange_b, damper.flange_b, inertia2.flange_a)]
@named model = ODESystem(connections,
t,
systems = [
fixed, inertia1, inertia2, spring, damper, speed_sensor,
- rel_speed_sensor, torque_sensor,
+ rel_speed_sensor, torque_sensor
])
sys = structural_simplify(model)
@@ -261,7 +259,7 @@ end
@test all(sol[rel_speed_sensor.w_rel.u] .== sol[speed_sensor.w.u])
@test all(sol[torque_sensor.tau.u] .== -sol[inertia1.flange_b.tau])
- prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, Pair[], (0, 10.0))
+ prob = DAEProblem(sys, D.(unknowns(sys)) .=> 0.0, Pair[], (0, 10.0))
sol = solve(prob, DFBDF())
@test SciMLBase.successful_retcode(sol)
@test all(sol[inertia1.w] .== 0)
diff --git a/test/Mechanical/translational.jl b/test/Mechanical/translational.jl
index 2af419ea3..5d4757b32 100644
--- a/test/Mechanical/translational.jl
+++ b/test/Mechanical/translational.jl
@@ -1,23 +1,23 @@
using ModelingToolkit, OrdinaryDiffEq, Test
+using ModelingToolkit: t, D
using ModelingToolkitStandardLibrary.Blocks
+import ModelingToolkitStandardLibrary: Mechanical
import ModelingToolkitStandardLibrary.Mechanical.Translational as TV
import ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition as TP
-
-@parameters t
-D = Differential(t)
+using DynamicQuantities: @u_str
@testset "Free" begin
function System(; name)
systems = @named begin
- acc = TV.Acceleration()
- a = Constant(; k = -10)
+ acc = TV.Acceleration(false)
+ a = Constant(; k = -10, output.unit = u"m/s^2")
mass = TV.Mass(; m = 100)
free = TV.Free()
end
eqs = [connect(a.output, acc.a)
- connect(mass.flange, acc.flange, free.flange)]
+ connect(mass.flange, acc.flange, free.flange)]
ODESystem(eqs, t, [], []; name, systems)
end
@@ -47,7 +47,7 @@ end
function simplify_and_solve(damping, spring, body, ground)
eqs = [connect(spring.flange_a, body.flange, damping.flange_a)
- connect(spring.flange_b, damping.flange_b, ground.flange)]
+ connect(spring.flange_b, damping.flange_b, ground.flange)]
@named model = ODESystem(eqs, t; systems = [ground, body, spring, damping])
@@ -83,15 +83,15 @@ end
@named gp = TP.Fixed(s_0 = 1)
@named fv = TV.Force()
- @named fp = TP.Force(use_support = false)
+ @named fp = TP.Force()
- @named source = Sine(frequency = 3, amplitude = 2)
+ @named source = Sine(frequency = 3, amplitude = 2, unit = u"N")
function System(damping, spring, body, ground, f, source)
eqs = [connect(f.f, source.output)
- connect(f.flange, body.flange)
- connect(spring.flange_a, body.flange, damping.flange_a)
- connect(spring.flange_b, damping.flange_b, ground.flange)]
+ connect(f.flange, body.flange)
+ connect(spring.flange_a, body.flange, damping.flange_a)
+ connect(spring.flange_b, damping.flange_b, ground.flange)]
@named model = ODESystem(eqs, t;
systems = [ground, body, spring, damping, f, source])
@@ -126,19 +126,19 @@ end
spring = TV.Spring(; k = 1000)
- src1 = Sine(frequency = 100, amplitude = 2)
- src2 = Sine(frequency = 100, amplitude = -1)
+ src1 = Sine(frequency = 100, amplitude = 2, output__unit = u"m")
+ src2 = Sine(frequency = 100, amplitude = -1, output__unit = u"N")
- pos_value = RealInput()
- force_output = RealOutput()
+ pos_value = RealInput(unit = u"m")
+ force_output = RealOutput(unit = u"N")
end
eqs = [connect(pos.s, src1.output)
- connect(force.f, src2.output)
- connect(spring.flange_a, pos.flange, force_sensor.flange)
- connect(spring.flange_b, force.flange, pos_sensor.flange)
- connect(pos_value, pos_sensor.output)
- connect(force_output, force_sensor.output)]
+ connect(force.f, src2.output)
+ connect(spring.flange_a, pos.flange, force_sensor.flange)
+ connect(spring.flange_b, force.flange, pos_sensor.flange)
+ connect(pos_value, pos_sensor.output)
+ connect(force_output, force_sensor.output)]
ODESystem(eqs, t, [], []; name, systems)
end
diff --git a/test/Mechanical/translational_modelica.jl b/test/Mechanical/translational_modelica.jl
index 015b36572..e5b17e923 100644
--- a/test/Mechanical/translational_modelica.jl
+++ b/test/Mechanical/translational_modelica.jl
@@ -1,19 +1,18 @@
using ModelingToolkit, OrdinaryDiffEq, Test
+using ModelingToolkit: t, D
using ModelingToolkitStandardLibrary.Blocks
-import ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica as TP
-
-@parameters t
-D = Differential(t)
+import ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica as TM
+using DynamicQuantities: @u_str
@testset "spring damper mass fixed" begin
- @named damper = TP.Damper(; d = 1)
- @named spring = TP.Spring(; c = 1, s_rel0 = 1)
- @named mass = TP.Mass(; m = 1, v = 1)
- @named fixed = TP.Fixed(s0 = 1)
+ @named damper = TM.Damper(; d = 1)
+ @named spring = TM.Spring(; c = 1, s_rel0 = 1)
+ @named mass = TM.Mass(; m = 1, v = 1)
+ @named fixed = TM.Fixed(s0 = 1)
eqs = [connect(spring.flange_a, mass.flange_a, damper.flange_a)
- connect(spring.flange_b, damper.flange_b, fixed.flange)]
+ connect(spring.flange_b, damper.flange_b, fixed.flange)]
@named model = ODESystem(eqs, t; systems = [fixed, mass, spring, damper])
@@ -29,17 +28,17 @@ D = Differential(t)
end
@testset "driven spring damper mass" begin
- @named damper = TP.Damper(; d = 1)
- @named spring = TP.Spring(; c = 1, s_rel0 = 1)
- @named mass = TP.Mass(; m = 1, v = 1)
- @named fixed = TP.Fixed(; s0 = 1)
- @named force = TP.Force(use_support = false)
- @named source = Sine(frequency = 3, amplitude = 2)
+ @named damper = TM.Damper(; d = 1)
+ @named spring = TM.Spring(; c = 1, s_rel0 = 1)
+ @named mass = TM.Mass(; m = 1, v = 1)
+ @named fixed = TM.Fixed(; s0 = 1)
+ @named force = TM.Force(use_support = false)
+ @named source = Sine(frequency = 3, amplitude = 2, output__unit = u"N")
eqs = [connect(force.f, source.output)
- connect(force.flange, mass.flange_a)
- connect(spring.flange_a, mass.flange_b, damper.flange_a)
- connect(spring.flange_b, damper.flange_b, fixed.flange)]
+ connect(force.flange, mass.flange_a)
+ connect(spring.flange_a, mass.flange_b, damper.flange_a)
+ connect(spring.flange_b, damper.flange_b, fixed.flange)]
@named model = ODESystem(eqs, t;
systems = [fixed, mass, spring, damper, force, source])
diff --git a/test/Thermal/demo.jl b/test/Thermal/demo.jl
index 0cc950c56..a03aee0d1 100644
--- a/test/Thermal/demo.jl
+++ b/test/Thermal/demo.jl
@@ -1,8 +1,7 @@
using ModelingToolkitStandardLibrary.Thermal, ModelingToolkit, OrdinaryDiffEq, Test
+using ModelingToolkit: t, D
using OrdinaryDiffEq: ReturnCode.Success
-
-@parameters t
-D = Differential(t)
+using DynamicQuantities: @u_str
# Modelica example
@testset "demo" begin
@@ -16,7 +15,7 @@ D = Differential(t)
connect(mass1.port, conduction.port_a),
connect(conduction.port_b, mass2.port),
connect(mass1.port, Tsensor1.port),
- connect(mass2.port, Tsensor2.port),
+ connect(mass2.port, Tsensor2.port)
]
@named model = ODESystem(connections, t,
diff --git a/test/Thermal/thermal.jl b/test/Thermal/thermal.jl
index d8f69fd41..1508a7646 100644
--- a/test/Thermal/thermal.jl
+++ b/test/Thermal/thermal.jl
@@ -1,9 +1,9 @@
using ModelingToolkitStandardLibrary.Thermal, ModelingToolkit, OrdinaryDiffEq, Test
+using ModelingToolkit: t, D
using ModelingToolkitStandardLibrary.Blocks: Constant, Step
using OrdinaryDiffEq: ReturnCode.Success
+using DynamicQuantities: @u_str
-@parameters t
-D = Differential(t)
#=
# Test HeatCapacitor, TemperatureSensor, RelativeTemperatureSensor, FixedTemperature
@testset "Heat systems" begin
@@ -210,21 +210,21 @@ end
@named T_core = TemperatureSensor()
@named convection = ConvectiveConductor(G = 25)
@named environment = PrescribedTemperature()
- @named amb = Constant(k = T_amb)
- @named core_losses_const = Constant(k = 500)
+ @named amb = Constant(k = T_amb, output.unit = u"K")
+ @named core_losses_const = Constant(k = 500, output.unit = u"")
@named winding_losses = Step(height = 900, offset = 100, start_time = 360,
duration = Inf, smooth = false)
connections = [connect(windingLosses.port, winding.port)
- connect(coreLosses.port, core.port)
- connect(winding.port, winding2core.port_a)
- connect(winding2core.port_b, core.port)
- connect(winding.port, T_winding.port)
- connect(core.port, T_core.port)
- connect(winding2core.port_b, convection.solid)
- connect(convection.fluid, environment.port)
- connect(amb.output, environment.T)
- connect(winding_losses.output, windingLosses.Q_flow)
- connect(core_losses_const.output, coreLosses.Q_flow)]
+ connect(coreLosses.port, core.port)
+ connect(winding.port, winding2core.port_a)
+ connect(winding2core.port_b, core.port)
+ connect(winding.port, T_winding.port)
+ connect(core.port, T_core.port)
+ connect(winding2core.port_b, convection.solid)
+ connect(convection.fluid, environment.port)
+ connect(amb.output, environment.T)
+ connect(winding_losses.output, windingLosses.Q_flow)
+ connect(core_losses_const.output, coreLosses.Q_flow)]
@named model = ODESystem(connections, t,
systems = [
@@ -240,7 +240,7 @@ end
@test sol.retcode == Success
@test sol[T_winding.T] == sol[winding.T]
@test sol[T_core.T] == sol[core.T]
- @test sol[-core.port.Q_flow] ==
+ @test sol[-core.port.Q_flow] ≈
sol[coreLosses.port.Q_flow + convection.solid.Q_flow + winding2core.port_b.Q_flow]
@test sol[T_winding.T][end] >= 500 # not good but better than nothing
@test sol[T_core.T] <= sol[T_winding.T]
diff --git a/test/chua_circuit.jl b/test/chua_circuit.jl
index c12896d3c..3cb8fdd5a 100644
--- a/test/chua_circuit.jl
+++ b/test/chua_circuit.jl
@@ -1,25 +1,27 @@
using ModelingToolkit
+using ModelingToolkit: t
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Electrical: OnePort
using OrdinaryDiffEq
using OrdinaryDiffEq: ReturnCode.Success
using IfElse: ifelse
+using DynamicQuantities: @u_str
@testset "Chua Circuit" begin
- @parameters t
-
- @component function NonlinearResistor(; name, Ga, Gb, Ve)
- @named oneport = OnePort()
- @unpack v, i = oneport
- pars = @parameters Ga=Ga Gb=Gb Ve=Ve
- eqs = [
+ @mtkmodel NonlinearResistor begin
+ @extend OnePort()
+ @parameters begin
+ Ga, [description = "Conductance in inner voltage range", unit = u"1/Ω"]
+ Gb, [description = "Conductance in outer voltage range", unit = u"1/Ω"]
+ Ve, [description = "Inner voltage range limit", unit = u"V"]
+ end
+ @equations begin
i ~ ifelse(v < -Ve,
Gb * (v + Ve) - Ga * Ve,
ifelse(v > Ve,
Gb * (v - Ve) + Ga * Ve,
- Ga * v)),
- ]
- extend(ODESystem(eqs, t, [], pars; name = name), oneport)
+ Ga * v))
+ end
end
@named L = Inductor(L = 18, i = 0.0)
@@ -33,14 +35,14 @@ using IfElse: ifelse
@named Gnd = Ground()
connections = [connect(L.p, G.p)
- connect(G.n, Nr.p)
- connect(Nr.n, Gnd.g)
- connect(C1.p, G.n)
- connect(L.n, Ro.p)
- connect(G.p, C2.p)
- connect(C1.n, Gnd.g)
- connect(C2.n, Gnd.g)
- connect(Ro.n, Gnd.g)]
+ connect(G.n, Nr.p)
+ connect(Nr.n, Gnd.g)
+ connect(C1.p, G.n)
+ connect(L.n, Ro.p)
+ connect(G.p, C2.p)
+ connect(C1.n, Gnd.g)
+ connect(C2.n, Gnd.g)
+ connect(Ro.n, Gnd.g)]
@named model = ODESystem(connections, t, systems = [L, Ro, G, C1, C2, Nr, Gnd])
sys = structural_simplify(model)
diff --git a/test/multi_domain.jl b/test/multi_domain.jl
index 6b113c698..501d75d46 100644
--- a/test/multi_domain.jl
+++ b/test/multi_domain.jl
@@ -5,12 +5,11 @@ using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkitStandardLibrary.Thermal
import ModelingToolkitStandardLibrary
using ModelingToolkit, OrdinaryDiffEq, Test
+using ModelingToolkit: t, D
using OrdinaryDiffEq: ReturnCode.Success
+using DynamicQuantities: @u_str
# using Plots
-@parameters t
-D = Differential(t)
-
@testset "DC motor" begin
R = 0.5
L = 4.5e-3
@@ -21,25 +20,25 @@ D = Differential(t)
tau_L_step = -3
@named ground = Ground()
@named source = Voltage()
- @named voltage_step = Blocks.Step(height = V_step, start_time = 0)
+ @named voltage_step = Blocks.Step(height = V_step, start_time = 0, output__unit = u"V")
@named R1 = Resistor(R = R)
@named L1 = Inductor(L = L, i = 0.0)
@named emf = EMF(k = k)
@named fixed = Fixed()
@named load = Torque(use_support = false)
- @named load_step = Blocks.Step(height = tau_L_step, start_time = 3)
+ @named load_step = Blocks.Step(height = tau_L_step, start_time = 3, output__unit = u"N*m")
@named inertia = Inertia(J = J)
@named friction = Damper(d = f)
connections = [connect(fixed.flange, emf.support, friction.flange_b)
- connect(emf.flange, friction.flange_a, inertia.flange_a)
- connect(inertia.flange_b, load.flange)
- connect(load_step.output, load.tau)
- connect(voltage_step.output, source.V)
- connect(source.p, R1.p)
- connect(R1.n, L1.p)
- connect(L1.n, emf.p)
- connect(emf.n, source.n, ground.g)]
+ connect(emf.flange, friction.flange_a, inertia.flange_a)
+ connect(inertia.flange_b, load.flange)
+ connect(load_step.output, load.tau)
+ connect(voltage_step.output, source.V)
+ connect(source.p, R1.p)
+ connect(R1.n, L1.p)
+ connect(L1.n, emf.p)
+ connect(emf.n, source.n, ground.g)]
@named model = ODESystem(connections, t,
systems = [
@@ -53,7 +52,7 @@ D = Differential(t)
load,
load_step,
inertia,
- friction,
+ friction
])
sys = structural_simplify(model)
@@ -72,7 +71,7 @@ D = Differential(t)
@test sol[inertia.w][idx_t]≈(dc_gain * [V_step; -tau_L_step])[2] rtol=1e-3
@test sol[emf.i][idx_t]≈(dc_gain * [V_step; -tau_L_step])[1] rtol=1e-3
- prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, Pair[], (0, 6.0))
+ prob = DAEProblem(sys, D.(unknowns(sys)) .=> 0.0, Pair[], (0, 6.0))
sol = solve(prob, DFBDF())
@test sol.retcode == Success
# EMF equations
@@ -103,27 +102,27 @@ end
tau_L_step = -3
@named ground = Ground()
@named source = Voltage()
- @named voltage_step = Blocks.Step(height = V_step, start_time = 0)
+ @named voltage_step = Blocks.Step(height = V_step, start_time = 0, output__unit = u"V")
@named R1 = Resistor(R = R)
@named L1 = Inductor(L = L, i = 0.0)
@named emf = EMF(k = k)
@named fixed = Fixed()
@named load = Torque(use_support = false)
- @named load_step = Blocks.Step(height = tau_L_step, start_time = 3)
+ @named load_step = Blocks.Step(height = tau_L_step, start_time = 3, output__unit = u"N*m")
@named inertia = Inertia(J = J)
@named friction = Damper(d = f)
@named speed_sensor = SpeedSensor()
connections = [connect(fixed.flange, emf.support, friction.flange_b)
- connect(emf.flange, friction.flange_a, inertia.flange_a)
- connect(inertia.flange_b, load.flange)
- connect(inertia.flange_b, speed_sensor.flange)
- connect(load_step.output, load.tau)
- connect(voltage_step.output, source.V)
- connect(source.p, R1.p)
- connect(R1.n, L1.p)
- connect(L1.n, emf.p)
- connect(emf.n, source.n, ground.g)]
+ connect(emf.flange, friction.flange_a, inertia.flange_a)
+ connect(inertia.flange_b, load.flange)
+ connect(inertia.flange_b, speed_sensor.flange)
+ connect(load_step.output, load.tau)
+ connect(voltage_step.output, source.V)
+ connect(source.p, R1.p)
+ connect(R1.n, L1.p)
+ connect(L1.n, emf.p)
+ connect(emf.n, source.n, ground.g)]
@named model = ODESystem(connections, t,
systems = [
@@ -138,7 +137,7 @@ end
load_step,
inertia,
friction,
- speed_sensor,
+ speed_sensor
])
sys = structural_simplify(model)
@@ -161,7 +160,7 @@ end
@test all(sol[inertia.w] .== sol[speed_sensor.w.u])
- prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, Pair[], (0, 6.0))
+ prob = DAEProblem(sys, D.(unknowns(sys)) .=> 0.0, Pair[], (0, 6.0))
sol = solve(prob, DFBDF())
@test sol.retcode == Success
@@ -183,15 +182,15 @@ end
@testset "El. Heating Circuit" begin
@named ground = Ground()
@named source = Voltage()
- @named voltage_sine = Blocks.Sine(amplitude = 220, frequency = 1)
+ @named voltage_sine = Blocks.Sine(amplitude = 220, frequency = 1, output__unit = u"V")
@named heating_resistor = HeatingResistor(R_ref = 100, alpha = 1e-3, T_ref = 293.15)
@named thermal_conductor = ThermalConductor(G = 50)
@named env = FixedTemperature(T = 273.15 + 20)
connections = [connect(source.n, ground.g, heating_resistor.n)
- connect(source.p, heating_resistor.p)
- connect(voltage_sine.output, source.V)
- connect(heating_resistor.heat_port, thermal_conductor.port_a)
- connect(thermal_conductor.port_b, env.port)]
+ connect(source.p, heating_resistor.p)
+ connect(voltage_sine.output, source.V)
+ connect(heating_resistor.heat_port, thermal_conductor.port_a)
+ connect(thermal_conductor.port_b, env.port)]
@named model = ODESystem(connections, t,
systems = [
@@ -200,12 +199,13 @@ end
source,
heating_resistor,
thermal_conductor,
- env,
+ env
])
sys = structural_simplify(model)
prob = ODEProblem(sys, Pair[], (0, 6.0))
sol = solve(prob, Rodas4())
@test sol.retcode == Success
+ ### TODO
@test sol[source.v * source.i] == -sol[env.port.Q_flow]
end
diff --git a/test/qa.jl b/test/qa.jl
new file mode 100644
index 000000000..5d580adc1
--- /dev/null
+++ b/test/qa.jl
@@ -0,0 +1,11 @@
+using ModelingToolkitStandardLibrary, Aqua
+@testset "Aqua" begin
+ Aqua.find_persistent_tasks_deps(ModelingToolkitStandardLibrary)
+ Aqua.test_ambiguities(ModelingToolkitStandardLibrary, recursive = false)
+ Aqua.test_deps_compat(ModelingToolkitStandardLibrary)
+ Aqua.test_piracies(ModelingToolkitStandardLibrary)
+ Aqua.test_project_extras(ModelingToolkitStandardLibrary)
+ Aqua.test_stale_deps(ModelingToolkitStandardLibrary)
+ Aqua.test_unbound_args(ModelingToolkitStandardLibrary)
+ Aqua.test_undefined_exports(ModelingToolkitStandardLibrary)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 178984894..282175ab7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,9 @@
using SafeTestsets
+@safetestset "Quality Assurance" begin
+ include("qa.jl")
+end
+
# Blocks
@safetestset "Blocks: math" begin
include("Blocks/math.jl")