add Unitful example to docs

docs/src/examples.md

@@ -29,6 +29,46 @@ using MonteCarloMeasurements, Unitful
 (1..2)u"m"
 ```
 
+### Example: Solar collector energy transfer
+The following example estimates the amount of thermal power transferred from a solar collector embedded in a concrete floor, to a water reservoir. The power is computed by measuring the temperature difference, $\Delta T$, between the solar collectors circulating warm water going into the collector tank and the colder returning water. Using the mass-flow rate and the specific heat capacity of water, we can estimate the power transfer. No flow meter is installed, so the flow is estimated and subject to large uncertainty.
+```@example solar
+using MonteCarloMeasurements
+using Unitful
+using Unitful: W, kW, m, mm, hr, K, g, J, l, s
+
+ΔT = (3.5 ± 0.8)K # The temperature difference varies slightly between different flow circuits.
+specific_heat_water = 4.19J/(g*K)
+density_water = 1e6g/m^3
+flow = 8*(1..2.5)*l/(60s) # 8 solar collector circuits, each with an estimated flow rate of between 1 and 2.5 l/minute
+mass_flow = flow * density_water |> upreferred # Water flow in mass per second
+power = uconvert(W, mass_flow * specific_heat_water * ΔT) # Power in Watts
+```
+
+Some power is lost to the ground in which the heat-exchanger circuits are embedded, we estimate this to be between 10 and 50% of the total power.
+```@example solar
+ground_losses = (0.1..0.5) * power # Between 10-50% power loss to ground
+reservoir_volume = 7m*3m*1.5m
+```
+
+The energy transfered during 6hrs solar collection can be estimated as
+```@example solar
+energy_6_hrs = (power - ground_losses)*6hr
+```
+and this energy transfer will increase the temperature in the reservoir by
+```@example solar
+ΔT_reservoir_6hr = energy_6_hrs/(reservoir_volume*density_water*specific_heat_water) |> upreferred
+```
+
+Finally, we visualize the distributions associated with the estimated quantities:
+```@example solar
+using Plots
+figh = plot(ΔT_reservoir_6hr, xlabel="\$ΔT [K]\$", ylabel="\$P(ΔT)\$", title="Temperature increase 6hrs sun")
+qs   = 0:0.01:1
+Qs   = pquantile.(ΔT_reservoir_6hr, qs)
+figq = plot(qs, Qs, xlabel="∫\$P(ΔT)\$")
+plot(figh, figq)
+```
+
 
 ## Monte-Carlo sampling properties
 The variance introduced by Monte-Carlo sampling has some fortunate and some unfortunate properties. It decreases as 1/N, where N is the number of particles/samples. This unfortunately means that to get half the standard deviation in your estimate, you need to quadruple the number of particles. On the other hand, this variance does not depend on the dimension of the space, which is very fortunate.

test/test_deconstruct.jl

@@ -5,12 +5,14 @@ using MonteCarloMeasurements: nakedtypeof, build_container, build_mutable_contai
 using ControlSystemsBase, Test, GenericSchur
 ControlSystemsBase.TransferFunction(matrix::Array{<:ControlSystemsBase.SisoRational,2}, Ts, ::Int64, ::Int64) = TransferFunction(matrix,Ts)
 
+Continuous = ControlSystemsBase.Continuous
+
 
 @testset "deconstruct" begin
     @info "Testing deconstruct"
     unsafe_comparisons()
     N = 50
-    P = tf(1 +0.1StaticParticles(N), [1, 1+0.1StaticParticles(N)])
+    P = ss(tf(1 +0.1StaticParticles(N), [1, 1+0.1StaticParticles(N)]))
     f = x->c2d(x,0.1)
     w = Workspace(f,P)
     @time Pd = w(P)
@@ -34,22 +36,17 @@ ControlSystemsBase.TransferFunction(matrix::Array{<:ControlSystemsBase.SisoRatio
     p = 1 ± 0.1
     @test mean_object(p) == pmean(p)
     @test mean_object([p,p]) == pmean.([p,p])
-    @test mean_object(P) ≈ tf(tf(1,[1,1])) atol=1e-2
+    @test mean_object(P) ≈ ss(tf(tf(1,[1,1]))) atol=1e-2
 
 
 
-    @test nakedtypeof(P) == TransferFunction
-    @test nakedtypeof(typeof(P)) == TransferFunction
-    @test typeof(P) == TransferFunction{ControlSystemsBase.Continuous, ControlSystemsBase.SisoRational{StaticParticles{Float64,N}}}
+    @test nakedtypeof(P) == StateSpace
+    @test nakedtypeof(typeof(P)) == StateSpace
     P2 = build_container(P)
-    @test typeof(P2) == TransferFunction{ControlSystemsBase.Continuous, ControlSystemsBase.SisoRational{Float64}}
-    @test typeof(build_mutable_container(P)) == TransferFunction{ControlSystemsBase.Continuous, ControlSystemsBase.SisoRational{Particles{Float64,N}}}
+    @test typeof(P2) == StateSpace{Continuous, Float64}
+    @test typeof(build_mutable_container(P)) == StateSpace{Continuous, Particles{Float64, N}}
     @test has_particles(P)
-    @test has_particles(P.matrix)
-    @test has_particles(P.matrix[1])
-    @test has_particles(P.matrix[1].num)
-    @test has_particles(P.matrix[1].num.coeffs)
-    @test has_particles(P.matrix[1].num.coeffs[1])
+    @test has_particles(P.A)
 
     # P = tf(1 ± 0.1, [1, 1±0.1])
     # @benchmark foreach(i->c2d($(tf(1.,[1., 1])),0.1), 1:N) # 1.7 ms 1.2 Mb
@@ -76,31 +73,31 @@ ControlSystemsBase.TransferFunction(matrix::Array{<:ControlSystemsBase.SisoRatio
         resultsetter = MonteCarloMeasurements.get_result_setter(Pres)
         @test all(1:paths[1][3]) do i
             buffersetter(P,P2,i)
-            P.matrix[1].num.coeffs[1].particles[i] == P2.matrix[1].num.coeffs[1] &&
-            P.matrix[1].den.coeffs[2].particles[i] == P2.matrix[1].den.coeffs[2]
+
+            P.A[1].particles[i] == P2.A[1] && P.A[1].particles[i] == P2.A[1]
             P2res = f(P2)
             resultsetter(Pres, P2res, i)
-            Pres.matrix[1].num.coeffs[1].particles[i] == P2res.matrix[1].num.coeffs[1] &&
-            Pres.matrix[1].den.coeffs[2].particles[i] == P2res.matrix[1].den.coeffs[2]
+
+            Pres.A[1].particles[i] == P2res.A[1] && Pres.A[1].particles[i] == P2res.A[1]
         end
     end
 
     @test mean_object(complex(1. ± 0.1, 1.)) isa ComplexF64
     @test mean_object(complex(1. ± 0.1, 1.)) ≈ complex(1,1) atol=1e-3
 
     Ps = MonteCarloMeasurements.make_scalar(P)
-    @test MonteCarloMeasurements.particletypetuple(Ps.matrix[1].num.coeffs[1]) == (Float64,1,Particles)
-    @test MonteCarloMeasurements.particletypetuple(MonteCarloMeasurements.restore_scalar(Ps,50).matrix[1].num.coeffs[1]) == (Float64,50,Particles)
+    @test MonteCarloMeasurements.particletypetuple(Ps.A[1]) == (Float64,1,Particles)
+    @test MonteCarloMeasurements.particletypetuple(MonteCarloMeasurements.restore_scalar(Ps,50).A[1]) == (Float64,50,Particles)
 
     unsafe_comparisons(false)
 
 
     N = 50
-    P = tf(1 +0.1StaticParticles(N), [1, 1+0.1StaticParticles(N)])
+    P = ss(tf(1 +0.1StaticParticles(N), [1, 1+0.1StaticParticles(N)]))
     f = x->c2d(x,0.1)
     res = MonteCarloMeasurements.array_of_structs(f, P)
     @test length(res) == N
-    @test res isa Vector{<:TransferFunction}
+    @test res isa Vector{<:StateSpace}
 
 
 end