Merge pull request #132 from acfr/hybridpassiveREN

Hybridpassive REN step 1: Incrementally Strictly Output Passive (ISOP) model

src/ParameterTypes/passive_ren.jl

@@ -9,11 +9,11 @@ mutable struct PassiveRENParams{T} <: AbstractRENParams{T}
     direct::DirectRENParams{T}
     αbar::T
     ν::T
-    # TODO: Add a field for incrementally strictly output passive model (ρ)
+    ρ::T
 end
 
 """
-    PassiveRENParams{T}(nu, nx, nv, ny, ν; <keyword arguments>) where T
+    PassiveRENParams{T}(nu, nx, nv, ny, ν, ρ; <keyword arguments>) where T
 
 Construct direct parameterisation of a passive REN.
 
@@ -22,8 +22,11 @@ Construct direct parameterisation of a passive REN.
 - `nx::Int`: Number of states.
 - `nv::Int`: Number of neurons.
 - `ny::Int`: Number of outputs.
-- `ν::Number=0`: Passivity parameter. Use ν>0 for incrementally strictly input passive model, and ν == 0 for incrementally passive model. 
-    
+- `ν::Number=0`: Passivity index. Use `ν > 0` for an incrementally strictly input passive model. Set both `ν = 0` and `ρ = 0` for incrementally passive model.
+- `ρ::Number=0`: Passivity index. Use `ρ > 0` for an incrementally strictly output passive model. 
+
+Note that setting both `ν,ρ > 0` or both `ν,ρ < 0` is not currently supported and will throw an error.
+
 # Keyword arguments
 
 - `nl::Function=relu`: Sector-bounded static nonlinearity.
@@ -35,7 +38,7 @@ See [`DirectRENParams`](@ref) for documentation of keyword arguments `init`, `ϵ
 See also [`GeneralRENParams`](@ref), [`ContractingRENParams`](@ref), [`LipschitzRENParams`](@ref).
 """
 function PassiveRENParams{T}(
-    nu::Int, nx::Int, nv::Int, ny::Int, ν::Number=T(0);
+    nu::Int, nx::Int, nv::Int, ny::Int, ν::Number=T(0), ρ::Number=T(0);
     nl::Function      = relu, 
     αbar::T           = T(1),
     init              = :random,
@@ -51,14 +54,23 @@ function PassiveRENParams{T}(
         error("Input and output must have the same dimension for passiveREN")
     end
 
+    # Check ρ and ν
+    if ρ*ν > 0
+        error("If ρ and ν are both positive, passiveREN could produce incorrect results. Please set at least one of them as zero. ")               
+    end
+
+    if ρ < 0 || ν < 0
+        @warn("Warning: negative passivity index detected, passivity is NOT guaranteed")
+    end 
+
     # Direct (implicit) params
     direct_ps = DirectRENParams{T}(
         nu, nx, nv, ny; 
         init, ϵ, bx_scale, bv_scale, polar_param, 
         D22_free=false, rng,
     )
 
-    return PassiveRENParams{T}(nl, nu, nx, nv, ny, direct_ps, αbar, ν)
+    return PassiveRENParams{T}(nl, nu, nx, nv, ny, direct_ps, αbar, ν, ρ)
 
 end
 
@@ -68,12 +80,12 @@ trainable(m::PassiveRENParams) = (direct = m.direct, )
 function direct_to_explicit(ps::PassiveRENParams{T}, return_h=false) where T
 
     # System sizes
-    nu = ps.nu
     ν = ps.ν
-
+    ρ = ps.ρ
+
     # Implicit parameters
     ϵ = ps.direct.ϵ
-    ρ = ps.direct.ρ
+    ρ_polar = ps.direct.ρ
     X = ps.direct.X
     polar_param = ps.direct.polar_param
 
@@ -87,28 +99,52 @@ function direct_to_explicit(ps::PassiveRENParams{T}, return_h=false) where T
     C2 = ps.direct.C2
     D21 = ps.direct.D21
 
-    # Constructing D22 for incrementally passive and incrementally strictly input passive. 
+    # Constructing D22 for incrementally (strictly input) passive and incrementally strictly output passive. 
     # See Eqns 31-33 of TAC paper 
-    # Currently converts to Hermitian to avoid numerical conditioning issues
     M = _M_pass(X3, Y3, ϵ)
 
-    D22 = ν*I + M
-    D21_imp = D21 - D12_imp'
+    if ρ == 0
+        # For ρ==0 case, I(SI)P model
+        D22 = ν*I + M
+        D21_imp = D21 - D12_imp'
+
+        𝑅  = _R_pass(D22, ν, ρ) 
+        Γ2 = _Γ2_pass(C2, D21_imp, B2_imp, 𝑅)
 
-    𝑅  = _R_pass(D22, ν)
-    Γ2 = _Γ2_pass(C2, D21_imp, B2_imp, 𝑅)
+        H = x_to_h(X, ϵ, polar_param, ρ_polar) + Γ2
+    else    
+        # For ρ!=0 case, ISOP model
+        D22 = _D22_pass(M, ρ) 
+        C2_imp = _C2_pass(D22, C2, ρ)
+        D21_imp = _D21_pass(D22, D21, D12_imp, ρ)
 
-    H = x_to_h(X, ϵ, polar_param, ρ) + Γ2
+        𝑅  = _R_pass(D22, ν, ρ)
 
+        Γ1 = _Γ1_pass(ps.nx, ps.ny, C2, D21, ρ, T) 
+        Γ2 = _Γ2_pass(C2_imp, D21_imp, B2_imp, 𝑅)
+
+        H = x_to_h(X, ϵ, polar_param, ρ_polar) + Γ2 - Γ1
+    end
+
     # Get explicit parameterisation
     !return_h && (return hmatrix_to_explicit(ps, H, D22))
     return H
 
 end
 
+_D22_pass(M, ρ) = ((I+M) \ I) / ρ   
+
+_C2_pass(D22, C2, ρ) = (D22'*(-2ρ*I) + I)*C2
+
+_D21_pass(D22, D21, D12_imp, ρ) = (D22'*(-2ρ*I) + I)*D21 - D12_imp'
+
 _M_pass(X3, Y3, ϵ) = X3'*X3 + Y3 - Y3' + ϵ*I
 
-_R_pass(D22, ν) = -2ν*I + D22 + D22'
+_R_pass(D22, ν, ρ) = -2ν*I + D22 + D22' + D22'*(-2ρ*I)*D22
+
+function _Γ1_pass(nx, ny, C2, D21, ρ, T) 
+    [C2'; D21'; zeros(T, nx, ny)] * (-2ρ*I) * [C2 D21 zeros(T, ny, nx)]
+end
 
 function _Γ2_pass(C2, D21_imp, B2_imp, 𝑅)
     [C2'; D21_imp'; B2_imp] * (𝑅 \ [C2 D21_imp B2_imp'])

test/ParameterTypes/passivity.jl

@@ -10,16 +10,17 @@ using Test
 rng = MersenneTwister(42)
 
 """
-Test passivity inequality
+Test incrementally strictly input passivity (isip)
 """
 batches = 100
 nu, nx, nv, ny = 6, 5, 10, 6
-ν= 1.0
+ν = 2.0
+ρ = 10.0
 T = 100
 
-# Test constructors
-ren_ps = PassiveRENParams{Float64}(nu, nx, nv, ny, ν; init=:random, rng)
-ren = REN(ren_ps)
+# Test constructor
+isip_ren_ps = PassiveRENParams{Float64}(nu, nx, nv, ny, ν, 0; init=:random, rng)
+isip_ren = REN(isip_ren_ps)
 
 # Different inputs with different initial conditions
 u0 = 10*randn(rng, nu, batches)
@@ -29,20 +30,43 @@ x0 = randn(rng, nx, batches)
 x1 = randn(rng, nx, batches)
 
 # Simulate
-x0n, y0 = ren(x0, u0)
-x1n, y1 = ren(x1, u1)
+x0n_isip, y0_isip = isip_ren(x0, u0)
+x1n_isip, y1_isip = isip_ren(x1, u1)
 
 # Dissipation condition
-ν = ren_ps.ν
-Q = zeros(ny, ny)
+ν = isip_ren_ps.ν
+Q_isip = zeros(ny, ny)
+R_isip = -2ν * Matrix(I, nu, nu)
 S = Matrix(I, nu, ny)
-R = -2ν * Matrix(I, nu, nu)
 
 # Test passivity
-M = [Q S'; S R]
-rhs = mat_norm2(M, vcat(y1 .- y0, u1 .- u0))
+M_isip = [Q_isip S'; S R_isip]
+P_isip = compute_p(isip_ren_ps)
+rhs_isip = mat_norm2(M_isip, vcat(y1_isip .- y0_isip, u1 .- u0))
+lhs_isip = mat_norm2(P_isip, x0n_isip .- x1n_isip) - mat_norm2(P_isip, x0 .- x1)
 
-P = compute_p(ren_ps)
-lhs = mat_norm2(P, x0n .- x1n) - mat_norm2(P, x0 .- x1)
+@test all(lhs_isip .<= rhs_isip)
 
-@test all(lhs .<= rhs)
+"""
+Test incrementally strictly output passivity (isop)
+"""
+isop_ren_ps = PassiveRENParams{Float64}(nu, nx, nv, ny, 0, ρ; init=:random, rng)
+isop_ren = REN(isop_ren_ps)
+
+# Simulate
+x0n_isop, y0_isop = isop_ren(x0, u0)
+x1n_isop, y1_isop = isop_ren(x1, u1)
+
+# Dissipation condition
+ρ = isop_ren_ps.ρ
+Q_isop = -2ρ * Matrix(I, ny, ny)
+R_isop = zeros(nu, nu)
+S = Matrix(I, nu, ny)
+
+# Test passivity
+M_isop = [Q_isop S'; S R_isop]
+P_isop = compute_p(isop_ren_ps)
+rhs_isop = mat_norm2(M_isop, vcat(y1_isop .- y0_isop, u1 .- u0))
+lhs_isop = mat_norm2(P_isop, x0n_isop .- x1n_isop) - mat_norm2(P_isop, x0 .- x1)
+
+@test all(lhs_isop .<= rhs_isop)