Merge pull request #274 from ven-k/vkb/guess

Rotational: default values are now guesses + fix tests

Author: ChrisRackauckas

.typos.toml

@@ -1,4 +1,5 @@
 [default.extend-words]
 Nd = "Nd"
 nin = "nin"
-coul = "coul"
+coul = "coul"
+isconnection = "isconnection"

docs/src/tutorials/dc_motor_pi.md

@@ -88,7 +88,7 @@ so that it can be represented as a system of `ODEs` (ordinary differential equat
 
 ```@example dc_motor_pi
 sys = structural_simplify(model)
-prob = ODEProblem(sys, [], (0, 6.0))
+prob = ODEProblem(sys, unknowns(sys) .=> 0.0, (0, 6.0))
 sol = solve(prob, Rodas4())
 
 p1 = Plots.plot(sol.t, sol[inertia.w], ylabel = "Angular Vel. in rad/s",
@@ -118,9 +118,11 @@ T(s) &= \dfrac{P(s)C(s)}{I + P(s)C(s)}
 
 ```@example dc_motor_pi
 using ControlSystemsBase
-matrices_S, simplified_sys = Blocks.get_sensitivity(model, :y)
+matrices_S, simplified_sys = Blocks.get_sensitivity(
+    model, :y, op = Dict(unknowns(sys) .=> 0.0))
 So = ss(matrices_S...) |> minreal # The output-sensitivity function as a StateSpace system
-matrices_T, simplified_sys = Blocks.get_comp_sensitivity(model, :y)
+matrices_T, simplified_sys = Blocks.get_comp_sensitivity(
+    model, :y, op = Dict(inertia.phi => 0.0, inertia.w => 0.0))
 To = ss(matrices_T...)# The output complementary sensitivity function as a StateSpace system
 bodeplot([So, To], label = ["S" "T"], plot_title = "Sensitivity functions",
     plotphase = false)
@@ -129,7 +131,8 @@ bodeplot([So, To], label = ["S" "T"], plot_title = "Sensitivity functions",
 Similarly, we may compute the loop-transfer function and plot its Nyquist curve
 
 ```@example dc_motor_pi
-matrices_L, simplified_sys = Blocks.get_looptransfer(model, :y)
+matrices_L, simplified_sys = Blocks.get_looptransfer(
+    model, :y, op = Dict(unknowns(sys) .=> 0.0))
 L = -ss(matrices_L...) # The loop-transfer function as a StateSpace system. The negative sign is to negate the built-in negative feedback
 Ms, ωMs = hinfnorm(So) # Compute the peak of the sensitivity function to draw a circle in the Nyquist plot
 nyquistplot(L, label = "\$L(s)\$", ylims = (-2.5, 0.5), xlims = (-1.2, 0.1),

src/Blocks/continuous.jl

@@ -12,6 +12,10 @@ Initial value of integrator state ``x`` can be set with `x`
 # Parameters:
 
   - `k`: Gain of integrator
+
+# Unknowns:
+
+  - `x`: State of Integrator. Defaults to 0.0.
 """
 @mtkmodel Integrator begin
     @extend u, y = siso = SISO()
@@ -26,6 +30,7 @@ Initial value of integrator state ``x`` can be set with `x`
         y ~ x
     end
 end
+
 """
     Derivative(; name, k = 1, T, x = 0.0)
 
@@ -49,6 +54,10 @@ A smaller `T` leads to a more ideal approximation of the derivative.
   - `k`: Gain
   - `T`: [s] Time constants (T>0 required; T=0 is ideal derivative block)
 
+# Unknowns:
+
+  - `x`: Unknown of Derivative. Defaults to 0.0.
+
 # Connectors:
 
   - `input`
@@ -160,8 +169,8 @@ Initial value of the state `x` can be set with `x`, and of derivative state `xd`
 @mtkmodel SecondOrder begin
     @extend u, y = siso = SISO()
     @variables begin
-        x(t) = 0.0, [description = "State of SecondOrder filter"]
-        xd(t) = 0.0, [description = "Derivative state of SecondOrder filter"]
+        x(t), [description = "State of SecondOrder filter", guess = 0.0]
+        xd(t), [description = "Derivative state of SecondOrder filter", guess = 0.0]
     end
     @parameters begin
         k = 1.0, [description = "Gain"]

src/Blocks/math.jl

@@ -99,7 +99,7 @@ Output difference between reference input (input1) and feedback input (input2).
 end
 
 """
-    Add(; name, k1 = 1, k2 = 1)
+    Add(; name, k1 = 1.0, k2 = 1.0)
 
 Output the sum of the two scalar inputs.
 
@@ -121,16 +121,16 @@ Output the sum of the two scalar inputs.
         output = RealOutput()
     end
     @parameters begin
-        k1 = 1, [description = "Gain of Add input1"]
-        k2 = 1, [description = "Gain of Add input2"]
+        k1 = 1.0, [description = "Gain of Add input1"]
+        k2 = 1.0, [description = "Gain of Add input2"]
     end
     @equations begin
         output.u ~ k1 * input1.u + k2 * input2.u
     end
 end
 
 """
-    Add(; name, k1 = 1, k2 = 1, k3 = 1)
+    Add(; name, k1 = 1.0, k2 = 1.0, k3 = 1.0)
 
 Output the sum of the three scalar inputs.
 
@@ -155,9 +155,9 @@ Output the sum of the three scalar inputs.
         output = RealOutput()
     end
     @parameters begin
-        k1 = 1, [description = "Gain of Add input1"]
-        k2 = 1, [description = "Gain of Add input2"]
-        k3 = 1, [description = "Gain of Add input3"]
+        k1 = 1.0, [description = "Gain of Add input1"]
+        k2 = 1.0, [description = "Gain of Add input2"]
+        k3 = 1.0, [description = "Gain of Add input3"]
     end
     @equations begin
         output.u ~ k1 * input1.u + k2 * input2.u + k3 * input3.u

src/Mechanical/Rotational/components.jl

@@ -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), [description = "Absolute rotation angle", guess = 0.0]
+        w(t), [description = "Absolute angular velocity", guess = 0.0]
+        a(t), [description = "Absolute angular acceleration", guess = 0.0]
     end
     @equations begin
         phi ~ flange_a.phi
@@ -86,7 +86,7 @@ Linear 1D rotational spring
 # Parameters:
 
   - `c`: [`N.m/rad`] Spring constant
-  - `phi_rel0`: [`rad`] Unstretched spring angle
+  - `phi_rel0`: [`rad`] Unstretched spring angle. Defaults to 0.0.
 """
 @mtkmodel Spring begin
     @extend phi_rel, tau = partial_comp = PartialCompliant()
@@ -156,7 +156,7 @@ Linear 1D rotational spring and damper
 
   - `d`: [`N.m.s/rad`] Damping constant
   - `c`: [`N.m/rad`] Spring constant
-  - `phi_rel0`: [`rad`] Unstretched spring angle
+  - `phi_rel0`: [`rad`] Unstretched spring angle. Defaults to 0.0
 """
 @mtkmodel SpringDamper begin
     @extend phi_rel, w_rel, tau = partial_comp = PartialCompliantWithRelativeStates()
@@ -206,8 +206,10 @@ This element characterizes any type of gear box which is fixed in the ground and
         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),
+        [description = "Relative angle between shaft a and the support", guess = 0.0]
+        phi_b(t),
+        [description = "Relative angle between shaft b and the support", guess = 0.0]
     end
     @equations begin
         phi_a ~ flange_a.phi - phi_support

src/Mechanical/Rotational/sensors.jl

@@ -81,7 +81,7 @@ Ideal sensor to measure the relative angular velocity
         w_rel = RealOutput()
     end
     @variables begin
-        phi_rel(t) = 0.0
+        phi_rel(t), [guess = 0.0]
     end
     @equations begin
         0 ~ flange_a.tau + flange_b.tau

src/Mechanical/Rotational/sources.jl

@@ -103,7 +103,7 @@ 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) [guess = 0.0] w(t) [guess = 0.0] a(t) [guess = 0.0]
     eqs = [phi ~ flange.phi - phi_support
            D(phi) ~ w]
     if exact

src/Mechanical/Rotational/utils.jl

@@ -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), [description = "Relative rotation angle between flanges", guess = 0.0]
+        tau(t), [description = "Torque between flanges", guess = 0.0]
     end
     @equations begin
         phi_rel ~ flange_b.phi - flange_a.phi
@@ -104,10 +104,11 @@ 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), [description = "Relative rotation angle between flanges", guess = 0.0]
+        w_rel(t), [description = "Relative angular velocity between flanges", guess = 0.0]
+        a_rel(t),
+        [description = "Relative angular acceleration between flanges", guess = 0.0]
+        tau(t), [description = "Torque between flanges", guess = 0.0]
     end
     @equations begin
         phi_rel ~ flange_b.phi - flange_a.phi
@@ -138,7 +139,8 @@ 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) [
+        description = "Absolute angle of support flange", guess = 0.0]
     if use_support
         @named support = Support()
         eqs = [support.phi ~ phi_support
@@ -173,7 +175,8 @@ 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) [
+        description = "Absolute angle of support flange", guess = 0.0]
     if use_support
         @named support = Support()
         eqs = [support.phi ~ phi_support

test/Blocks/continuous.jl

@@ -82,7 +82,7 @@ end
     @named pt2 = SecondOrder(; k = k, w = w, d = d)
     @named iosys = ODESystem(connect(c.output, pt2.input), t, systems = [pt2, c])
     sys = structural_simplify(iosys)
-    prob = ODEProblem(sys, Pair[], (0.0, 100.0))
+    prob = ODEProblem(sys, [unknowns(sys) .=> 0.0...; pt2.xd => 0.0], (0.0, 100.0))
     sol = solve(prob, Rodas4())
     @test sol.retcode == Success
     @test sol[pt2.output.u]≈pt2_func.(sol.t, k, w, d) atol=1e-3

test/Blocks/nonlinear.jl

@@ -39,7 +39,7 @@ using OrdinaryDiffEq: ReturnCode.Success
             systems = [source, lim, int])
         sys = structural_simplify(iosys)
 
-        prob = ODEProblem(sys, Pair[], (0.0, 10.0))
+        prob = ODEProblem(sys, unknowns(sys) .=> 0.0, (0.0, 10.0))
 
         sol = solve(prob, Rodas4())
         @test sol.retcode == Success