WIP: Power Converter Benchmark

.gitignore

@@ -19,4 +19,6 @@ paper/paper.pdf
 *.out
 *.bbl
 *.synctex.gz
-*.log
+*.log
+Control Systems/v-infinite-horizon-pwer-converter/sdp.py
+Control Systems/v-infinite-horizon-pwer-converter/vtailcost.py

Project.toml

@@ -5,6 +5,7 @@ version = "0.1.2"
 
 [deps]
 ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d"
+BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiscreteMarkovChains = "8abcb7ef-b365-4f7b-ac38-56893fb62f9f"
@@ -18,6 +19,7 @@ IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 IntervalLinearAlgebra = "92cbe1ac-9c24-436b-b0c9-5f7317aedcd5"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+LatexPrint = "d2208f48-c256-5759-9089-c25ed2a93924"
 LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
 Libtiff_jll = "89763e89-9b03-5906-acba-b20f662cd828"
 LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179"

docs/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
 CDDLib = "3391f64e-dcde-5f30-b752-e11513730f60"
 Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"

docs/src/examples/solvers/Power Converter.jl

@@ -0,0 +1,159 @@
+using Test     #src
+# # Example: DC-DC converter solved by [Naive abstraction](https://github.com/dionysos-dev/Dionysos.jl/blob/master/docs/src/manual/manual.md#solvers).
+#
+#md # [![Binder](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/DC-DC converter.ipynb)
+#md # [![nbviewer](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/DC-DC converter.ipynb)
+#
+# We consider a boost DC-DC converter which has been widely studied from the point of view of hybrid control, see for example in  [1, V.A],[2],[3].
+# This is a **safety problem** for a **switching system**.
+#
+# ![Boost DC-DC converter.](https://github.com/dionysos-dev/Dionysos.jl/blob/master/docs/assets/dcdcboost.jpg?raw=true)
+#
+# The state of the system is given by $x(t) = \begin{bmatrix} i_l(t) & v_c(t) \end{bmatrix}^\top$.
+# The switching system has two modes consisting in two-dimensional affine dynamics:
+# ```math
+# \dot{x} = f_p(x) = A_p x + b_p,\quad p=1,2
+# ```
+# with
+# ```math
+# A_1 = \begin{bmatrix} -\frac{r_l}{x_l} &0 \\ 0 & -\frac{1}{x_c}\frac{1}{r_0+r_c}  \end{bmatrix}, A_2= \begin{bmatrix} -\frac{1}{x_l}\left(r_l+\frac{r_0r_c}{r_0+r_c}\right) & -\frac{1}{x_l}\frac{r_0}{r_0+r_c}  \\ \frac{1}{x_c}\frac{r_0}{r_0+r_c}   & -\frac{1}{x_c}\frac{1}{r_0+r_c}  \end{bmatrix}, b_1 = b_2 = \begin{bmatrix} \frac{v_s}{x_l}\\0\end{bmatrix}.
+# ```
+# The goal is to design a controller to keep the state of the system in a safety region around the reference desired value, using as input only the switching
+# signal. In order to study the concrete system and its symbolic abstraction in a unified framework, we will solve the problem
+# for the sampled system with a sampling time $\tau$. For the construction of the relations in the abstraction, it is necessary to over-approximate attainable sets of
+# a particular cell. In this example, we consider the use of a growth bound function  [4, VIII.2, VIII.5] which is one of the possible methods to over-approximate
+# attainable sets of a particular cell based on the state reach by its center.
+#
+
+# First, let us import [StaticArrays](https://github.com/JuliaArrays/StaticArrays.jl) and [Plots](https://github.com/JuliaPlots/Plots.jl).
+
+import CDDLib
+import GLPK
+import OSQP
+using JuMP
+import Pavito
+import HiGHS
+import Ipopt
+
+# At this point, we import the useful Dionysos sub-modules.
+using Dionysos
+const DI = Dionysos
+const UT = DI.Utils
+const DO = DI.Domain
+const ST = DI.System
+const SY = DI.Symbolic
+const OP = DI.Optim
+const AB = OP.Abstraction
+
+# ### Definition of the system
+# we can import the module containing the DCDC problem like this 
+#include(joinpath(dirname(dirname(pathof(Dionysos))), "problems", "pc-osqp.jl"))
+include(joinpath(dirname(dirname(pathof(Dionysos))), "problems", "pc-osqp-rmaps.jl"))
+
+problem = PCOSQPRM.problem(CDDLib.Library(), Float64);
+
+# Finally, we select the method presented in [2] as our optimizer
+
+qp_solver = optimizer_with_attributes(
+    OSQP.Optimizer,
+    "eps_abs" => 1e-8,
+    "eps_rel" => 1e-8,
+    "max_iter" => 100000,
+    MOI.Silent() => true,
+);
+
+mip_solver = optimizer_with_attributes(HiGHS.Optimizer, MOI.Silent() => true);
+
+cont_solver = optimizer_with_attributes(Ipopt.Optimizer, MOI.Silent() => true);
+
+miqp_solver = optimizer_with_attributes(
+    Pavito.Optimizer,
+    "mip_solver" => mip_solver,
+    "cont_solver" => cont_solver,
+    MOI.Silent() => true,
+);
+
+algo = optimizer_with_attributes(
+    OP.BemporadMorari.Optimizer{Float64},
+    "continuous_solver" => qp_solver,
+    "mixed_integer_solver" => miqp_solver,
+    "indicator" => false,
+    "log_level" => 0,
+);
+
+# and use it to solve the given problem, with the help of the abstraction layer
+# MathOptInterface provided by [JuMP](https://github.com/jump-dev/JuMP.jl)
+optimizer = MOI.instantiate(algo)
+MOI.set(optimizer, MOI.RawOptimizerAttribute("problem"), problem)
+MOI.optimize!(optimizer)
+
+# We check the solver time
+MOI.get(optimizer, MOI.SolveTimeSec())
+
+# the termination status
+termination = MOI.get(optimizer, MOI.TerminationStatus())
+
+# the objective value
+objective_value = MOI.get(optimizer, MOI.ObjectiveValue())
+
+##@test objective_value ≈ 11.38 atol = 1e-2     #src
+
+# and recover the corresponding continuous trajectory
+xu = MOI.get(optimizer, ST.ContinuousTrajectoryAttribute());
+
+# ## A little bit of data visualization now:
+
+using Plots
+using Polyhedra
+using HybridSystems
+using Suppressor
+
+#Initialize our canvas
+fig = plot(;
+    aspect_ratio = :equal,
+    xtickfontsize = 10,
+    ytickfontsize = 10,
+    guidefontsize = 16,
+    titlefontsize = 14,
+);
+xlims!(-10.5, 3.0)
+ylims!(-10.5, 3.0)
+
+#Plot the discrete modes
+for mode in states(problem.system)
+    t =
+        (problem.system.ext[:modes] in [mode, mode + 11]) ? "XT" :
+        (
+            mode == problem.system.ext[:q_A] ? "A" :
+            (
+                mode == problem.system.ext[:q_B] ? "B" :
+                mode <= 11 ? string(mode) : string(mode - 11)
+            )
+        )
+    set = stateset(problem.system, mode)
+    plot!(set; color = :white)
+    UT.text_in_set_plot!(fig, set, t)
+end
+
+#Plot obstacles
+for i in eachindex(problem.system.ext[:obstacles])
+    set = problem.system.ext[:obstacles][i]
+    plot!(set; color = :black, opacity = 0.5)
+    UT.text_in_set_plot!(fig, set, "O$i")
+end
+
+#Plot trajectory
+x0 = problem.initial_set[2]
+x_traj = [x0, xu.x...]
+plot!(fig, UT.DrawTrajectory(x_traj));
+
+#Plot initial point
+plot!(fig, UT.DrawPoint(x0); color = :blue)
+annotate!(fig, x0[1], x0[2] - 0.5, "x0")
+
+# ### References
+#
+# 1. Gol, E. A., Lazar, M., & Belta, C. (2013). Language-guided controller synthesis for linear systems. IEEE Transactions on Automatic Control, 59(5), 1163-1176.
+# 1. Bemporad, A., & Morari, M. (1999). Control of systems integrating logic, dynamics, and constraints. Automatica, 35(3), 407-427.
+# 1. Legat B., Bouchat J., Jungers R. M. (2021). Abstraction-based branch and bound approach to Q-learning for hybrid optimal control. 3rd Annual Learning for Dynamics & Control Conference, 2021.
+# 1. Legat B., Bouchat J., Jungers R. M. (2021). Abstraction-based branch and bound approach to Q-learning for hybrid optimal control. https://www.codeocean.com/. https://doi.org/10.24433/CO.6650697.v1.

problems/pc-osqp-rmaps.jl

@@ -0,0 +1,203 @@
+module PCOSQPRM
+using Test
+
+using Dionysos: Dionysos
+const DI = Dionysos
+const UT = DI.Utils
+const PR = DI.Problem
+
+using FillArrays
+using Polyhedra
+using MathematicalSystems, HybridSystems
+using SemialgebraicSets
+using LinearAlgebra
+
+using CDDLib: CDDLib
+
+function system(lib, T::Type)
+
+    # Generate the automaton nodes
+    possible_values = [[-1, 0, 1], [-1, 0, 1], [-1, 0, 1], [0, 1], [0, 1], [0, 1]]
+    nvars = length(possible_values)
+
+    ## generate all possible combinations of the values
+    function generate_combinations(possible_values, depth)
+        if depth == 1
+            return [[v] for v in possible_values[1]]
+        end
+        return [
+            [v; c...] for v in possible_values[1] for c in generate_combinations(possible_values[2:end], depth - 1)
+        ]
+    end
+
+    ## generate all possible combinations of the values
+    nodes = generate_combinations(possible_values, nvars)
+    nnodes = length(nodes)
+
+    ## domains
+    function rect(x_l, x_u)
+        r = HalfSpace([-1], -T(x_l)) ∩ HalfSpace([1], T(x_u))
+        return polyhedron(r, lib)
+    end
+    function rect(x_l, x_u, y_l, y_u)
+        r =
+            HalfSpace([-1, 0], -T(x_l)) ∩ HalfSpace([1, 0], T(x_u)) ∩
+            HalfSpace([0, -1], -T(y_l)) ∩ HalfSpace([0, 1], T(y_u))
+        return polyhedron(r, lib)
+    end
+
+    pX = rect(-10,10)
+    pU = rect(-1, 1)
+
+	#Automaton:
+	automaton = GraphAutomaton(nnodes)
+    indexB = Int[]
+
+	#Resetmaps: (guards + reset)
+	A = T[
+		0.9999981424449595 9.91220860939534e-12 3.4481569821196496e-7 9.205526977551878e-5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+		-9.91220860939534e-12 0.9999981424449595 -9.205526977551878e-5 3.4481569821196496e-7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+		2.1535327064414276e-7 -2.644213073804883e-12 0.9999999080158245 -2.455696808386769e-5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+		2.644213073804883e-12 2.1535327064414276e-7 2.455696808386769e-5 0.9999999080158245 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+		0.0 0.0 0.0 0.0 0.9999999996930378 -2.477749999746475e-5 0.0 0.0 0.0 0.0 0.0 0.0
+		0.0 0.0 0.0 0.0 2.477749999746475e-5 0.9999999996930378 0.0 0.0 0.0 0.0 0.0 0.0
+		0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+		0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+		0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+		0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.99875 0.0 0.0
+		0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0012499999999999734 0.99875 0.0
+		0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
+	]
+
+	B = T[
+		6.256252307659102e-5 -3.128126153811649e-5 -3.128126153847453e-5 0.0 0.0 0.0
+		-2.0671108273679513e-16 5.418073430928138e-5 -5.418073430907468e-5 0.0 0.0 0.0
+		6.7365306173215874e-12 -3.368357184796113e-12 -3.3681734325413444e-12 0.0 0.0 0.0
+		1.0608933096199221e-16 5.833953603352774e-12 -5.8340596926953345e-12 0.0 0.0 0.0
+		0.0 0.0 0.0 0.0 0.0 0.0
+		0.0 0.0 0.0 0.0 0.0 0.0
+		1.0 0.0 0.0 0.0 0.0 0.0
+		0.0 1.0 0.0 0.0 0.0 0.0
+		0.0 0.0 1.0 0.0 0.0 0.0
+		0.0 0.0 0.0 0.013888888888888593 0.013888888888888593 0.013888888888888593
+		0.0 0.0 0.0 0.0 0.0 0.0
+		0.0 0.0 0.0 0.0 0.0 0.0
+	]
+
+    C = T[
+        1.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 
+        0.0 1.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 
+        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -5.5 5.5 
+    ]
+
+	#ResetMaps = Vector{ConstrainedLinearControlDiscreteSystem}(undef, length(domains)^2)
+	##cU = polyhedron(HalfSpace(T[0, 0, -1], 2) ∩ HalfSpace(T[0, 0, 1], 2), lib)
+
+	function is_enabled(from, to)
+        from_node = nodes[from]
+        to_node = nodes[to]
+
+        # On the first 3 elements of the nodes, we cannot go from -1 to 1 or 1 to -1
+        for i in 1:3
+            if from_node[i] == -1 && to_node[i] == 1
+                return false
+            end
+            if from_node[i] == 1 && to_node[i] == -1
+                return false
+            end
+        end
+
+        return true
+	end
+
+	ntransitions = 0
+	function maybe_add_transition(from, to)
+        #global ntransitions
+		if is_enabled(from, to)
+			ntransitions += 1
+			add_transition!(automaton, from, to, ntransitions)
+            append!(indexB, from)
+		end
+	end
+	for from in 1:nnodes, to in 1:nnodes
+		maybe_add_transition(from, to)
+	end
+
+    #modes = [ConstrainedContinuousIdentitySystem(12, i) for i in 1:nnodes]
+    modes = [ConstrainedContinuousIdentitySystem(12, FullSpace()) for i in 1:nnodes]
+
+    #L290: ERROR: DimensionMismatch: matrix A has dimensions (6,6), vector B has length 0
+    #resetmaps = [ConstrainedAffineMap(A, B*nodes[indexB[i]], FullSpace()) for i in 1:ntransitions]
+    #resetmaps = Fill(ConstrainedAffineMap(A, B*nodes[indexB[i]], FullSpace()), ntransitions)
+	resetmaps = Fill(ConstrainedLinearControlMap(A, B, FullSpace(), pU), ntransitions)
+
+	system = HybridSystem(
+		automaton,
+		# Modes
+		modes,
+		# Reset maps
+		resetmaps,
+        # Switching
+		Fill(ControlledSwitching(), ntransitions),
+        #Fill(AutonomousSwitching(), ntransitions),
+	)
+
+	system.ext[:modes] = nmodes(system)
+    system.ext[:q_C] = C
+    system.ext[:transitions] = ntransitions
+
+	return system
+end
+
+""""
+	problem(lib, T::Type, gamma=0.95, x_0 = [1.0, -6.0], N = 11, tail_cost::Bool = false)
+
+This function create the system with `PCOSQPRM.system`.
+
+Then, we define initial conditions (continuous and discrete states) to this system
+and set `N` as horizon, i.e., the number of allowed time steps.
+
+We instantiate our Optimal Control Problem by defining the state and tail costs.
+Notice that `state_cost` is defined to be `sum(gamma^t * norm(Cx_t))` for each time step `t`
+and the `tail_cost` is defined to be `gamma^N V(x_N)` that we considering constant in this version.
+
+Notice that we used `Fill` for all `N` time steps as we consider time-invariant costs.
+"""
+function problem(
+	lib = CDDLib.Library(),
+	T::Type = Float64;
+	gamma = 0.95,
+    q0 = 1,
+	x_0 = [1 for _ in 1:12], #The first state is the initial state #[0 for i in 1:12],
+	N = 3,
+	tail_cost::Bool = false,
+)
+	sys = system(lib, T)
+	#if tail_cost
+	#	transition_cost = Fill(UT.ZeroFunction(), nmodes(sys))
+    #end 
+
+
+	#state_cost = UT.QuadraticControlFunction(sys.ext[:q_C]' * sys.ext[:q_C])
+    state_cost = UT.QuadraticControlFunction(ones(T, 6, 6))
+
+    #x_0 = [0 for _ in 1:sys.ext[:modes]]
+    #problem = PR.OptimalControlProblem(sys, (q0, x_0), sys.ext[:modes], nothing, nothing, N)
+
+    #==#
+	problem = PR.OptimalControlProblem(
+		sys,
+		(q0, x_0),
+		sys.ext[:modes],
+        #state_cost,
+		Fill(Fill(state_cost, sys.ext[:modes]), N),#Fill(Fill(state_cost, sys.ext[:modes]),N),
+        #transition_cost,
+		Fill(Fill(UT.ZeroFunction(), sys.ext[:transitions]), N),
+		N,
+	)
+    #==#
+	return problem
+end
+
+end
+