{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "angry-indiana", "metadata": {}, "outputs": [], "source": [ "#[1]: Guillaume O. Berger, Raphaël M. Jungers, and Zheming Wang. Chance-constrained quasi-convex optimization with application to data-driven switched systems control.arXiv preprintarXiv:2101.01415, 2021.\n", "#[2]: Anne Rubbens. Data-driven stability of black-box switched linear systems. Master Thesis, 2021.\n", "\n", "#Packages required to use the notebook.\n", "using Random\n", "using LinearAlgebra\n", "using JuMP\n", "using CSDP\n", "using IncGammaBeta\n", "using HybridSystems\n", "using SwitchOnSafety\n", "using Plots\n", "using DelimitedFiles\n", "using NLsolve\n", "using Statistics\n", "using Combinatorics\n", "using Ipopt\n", "using DynamicPolynomials\n", "using PyPlot" ] }, { "cell_type": "code", "execution_count": 2, "id": "legislative-accuracy", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "generatetrajectories (generic function with 1 method)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# function generatesystem(m,n)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function generates a random system of m modes of dimension nxn and of unitary Euclidean norm \n", "#\n", "# Inputs: -m: Positive integer describing the number of modes.\n", "# -n: Positive integer describing the dimension of the modes.\n", "#\n", "# Output: -A: Matrix containing m random float matrices of dimension nxn and of unitary Euclidean norm.\n", "function generatesystem(m,n)\n", " A=[randn(n,n) for i=1:m]\n", " for i=1:m\n", " A[i]=A[i]/opnorm(A[i])\n", " end\n", " return A\n", "end\n", "# function generatetrajectories(m,n,l,A,N)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function generates N initial states x0_i drawn randomly, independently and uniformly on the sphere; \n", "# and N final states x of trajectories of length l associated to a system A and to the initial state x0_i: \n", "# x_i=A_(i_l)*A_(i_(l-1))*...*A_(i_1)*x0_i, where the modes A_(i_j) are drawn randomly, uniformly and\n", "# independently from A.\n", "#\n", "# Inputs: -m: Positive integer describing the number of modes of the system.\n", "# -n: Positive integer describing the dimension of the modes of the system.\n", "# -l: Positive integer describing the length of the trajectory generated.\n", "# -A: Matrix containing m float matrices of dimension nxn: the m modes A_1, ..., A_m of dimension nxn of the system.\n", "# -N: Positive integer describing the number of trajectories generated.\n", "#\n", "# Output: -x0: Float matrix of dimension nxN containing N initial states drawn randomly and uniformly on the sphere.\n", "# -x: Float matrix of dimension nxN containing the N final states of trajectories of length l associated to \n", "# the system A and to the initial states x0_i: x_i=A_(i_l)*A_(i_(l-1))*...*A_(i_1)*x0_i, where the modes\n", "# A_(i_j) are drawn randomly, uniformly and independently from A.\n", "function generatetrajectories(m,n,l,A,N)\n", "x0=zeros(n,N)\n", "x=zeros(n,N)\n", "for i=1:N\n", " y0=randn(n)\n", " x0[:,i]=y0/norm(y0)\n", " x[:,i]=x0[:,i]\n", " b=rand((1:m),l)\n", " for j=1:l\n", " x[:,i]=A[b[j]]*x[:,i]\n", " end\n", "end\n", "return x0, x\n", "end" ] }, { "cell_type": "code", "execution_count": 12, "id": "ready-saturn", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "trueJSR (generic function with 2 methods)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# function trueJSR(A,d)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function computes the JSR of system A using a Sum-Of-Square Lyapunov function of degree 2d.\n", "#\n", "# Inputs: -A: Float matrix of dimension mxnxn containing the m modes A_1, ..., A_m of dimension nxn of the system.\n", "# -d: Positive integer. 2d is the degree of the function used \n", "# Output: -jsr: Float containing the JSR of the system.\n", "\n", "function trueJSR(A,d=2)\n", " s = discreteswitchedsystem(A)\n", " optimizer_constructor = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)\n", " soslyapb(s, d, optimizer_constructor=optimizer_constructor, tol=1e-4, verbose=1)\n", " seq = sosbuildsequence(s, d, p_0=:Primal)\n", " psw = findsmp(seq)\n", " return psw.growthrate\n", "end " ] }, { "cell_type": "code", "execution_count": 4, "id": "aggregate-samuel", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "transfo (generic function with 1 method)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# function dlift(x,d,n)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function computes the d-lift of a vector x of length n.\n", "#\n", "# Inputs: -x: Float vector of dimension n.\n", "# -d: Positive integer describing the dimension of the d-lift.\n", "# -n: Positive integer describing the dimension of the vector x. \n", "#\n", "# Output: -X: Float vector of dimension D=binomial(n+d-1,d) containing the d-lift of the vector x.\n", "\n", "function dlift(x,d,n)\n", " D=binomial(n+d-1,d)\n", " X=ones(D)\n", " alpha=weakcompo(d,n)\n", " for i=1:D\n", " X[i]=sqrt(factorial(d))\n", " for j=1:n\n", " X[i]=X[i]*x[j]^(alpha[i,j])/sqrt(factorial(alpha[i,j]))\n", " end\n", " end\n", " return X\n", "end\n", "\n", "# function weakcompo(d,n)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function computes all the possible vectors alpha of length n, such that alpha1+...+alphan=d.\n", "#\n", "# Inputs: -d: Positive integer describing the norm 1 of the vector alpha.\n", "# -n: Positive integer describing the dimension of the vector alpha.\n", "#\n", "# Output: -aplha: Integer vector of dimension Dxn where D=binomial(n+d-1,d), containing all vectors alpha of length n\n", "# such that alpha1+...+alphan=d.\n", "function weakcompo(d,n)\n", " D=binomial(n+d-1,d)\n", " alpha=zeros(D,n)\n", " alpha1=zeros(D,n)\n", " for i=1:n\n", " alpha[i,i]=1\n", " end\n", " for k=2:d\n", " index=1\n", " for i=1:binomial(n+k-2,k-1)\n", " for j=1:n\n", " X=alpha[i,:]\n", " X[j]=X[j]+1\n", " c=true\n", " for l=1:index-1\n", " if X==alpha1[l,:]\n", " c=false\n", " end\n", " end\n", " if c==true\n", " alpha1[index,:]=X\n", " index+=1\n", " end\n", " end\n", " end\n", " alpha=alpha1\n", " alpha1=zeros(D,n)\n", " end\n", " return alpha\n", "end\n", "\n", "# function dliftmat(L,d,n)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function computes the d-lift of a matrix L of dimension n.\n", "#\n", "# Inputs: -L: Float matrix of dimension n.\n", "# -d: Positive integer describing the dimension of the d-lift.\n", "# -n: Positive integer describing the dimension of the vector x. \n", "#\n", "# Output: -Ld: Float matrix of dimension D=binomial(n+d-1,d) containing the d-lift of the matrix L.\n", "\n", "function dliftmat(L,d,n)\n", " A=zeros(1,n)\n", " B=zeros(1,n)\n", " D=binomial(n+d-1,d)\n", " alpha=weakcompo(d,n)\n", " M=zeros(d,d)\n", " Ld=ones(D,D)\n", " for i=1:D\n", " for j=1:D\n", " A=alpha[i,:]\n", " B=alpha[j,:]\n", " den=sqrt(prod(m -> factorial(A[m]), 1:n)*prod(m -> factorial(B[m]), 1:n))\n", " C=transfo(A,d,n)\n", " E=transfo(B,d,n)\n", " y=collect(permutations([i for i=1:d]))\n", " Ld[i,j]=sum(prod(L[convert(Int64,C[k]),convert(Int64,E[y[p][k]])] for k=1:d) for p=1:length(y))/den\n", " end\n", " end\n", " return Ld\n", "end\n", "function transfo(x,d,n)\n", " x1=x\n", " X=zeros(d,1)\n", " index=1\n", " i=1\n", " while index<=d\n", " while x1[i]>0\n", " X[index]=i\n", " x1[i]=x1[i]-1\n", " index+=1\n", " end\n", " i=i+1\n", " end\n", " return X\n", "end" ] }, { "cell_type": "code", "execution_count": 16, "id": "gorgeous-animation", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "normmax (generic function with 1 method)" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# function proba_bounds(m,n,l,N,beta,beta1,tol,d,x0,x,C)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function computes probabilistic upper bounds on the JSR of a system of m modes of dimension nxn, based\n", "# on N trajectories of length l of initial state x0_i and final state x_i for i=1...N. The upper bounds hold \n", "# with confidence at least beta1+beta-1.\n", "# The function computes an approximate CQLF or SOS Lyapunov function of degree 2d by solving the problem \n", "# min gamma s.t. (x_i)^[d]' P (x_i)^[d] ≤ gamma^(2dl) (x0_i)^[d]' P (x0_i)^[d] for i=1...N; P≥I; ||P||≤ C (1)\n", "# where x^[d] is the d-lift of x. The problem is solved by bisection over gamma, with an accuracy tol.\n", "# If d=2, the fonction computes the probabilistic bounds given in [1, Cor. 6] (based on a chance-constrained\n", "# approach), [2, Thm. 3.1] and [2,Thm. 3.2] (based on sensivity analysis). Else, the function computes the upper bound given in \n", "# [2, Thm. 3.3] (based on sensivity analysis).\n", "# For the bounds of [2], two sampled problems have to be solved: (1) and the computation of the maximal norm of the \n", "# modes of the system. beta_1 is the confidence level of the probabilistic upper bound on the maximal norm of the \n", "# modes, while beta is the confidence level on the probabilistic bound on the JSR of the system based on Problem (1)\n", "# and involving the maximal norm of the modes, described in [2, Prop. 3]. For all bounds, the total confidence\n", "# level is fixed to beta+beta1-1.\n", "#\n", "# Inputs: -m: Positive integer describing the number of modes of the system.\n", "# -n: Positive integer describing the dimension of the modes of the system.\n", "# -l: Positive integer describing the length of the trajectories.\n", "# -N: Positive integer describing the number of trajectories sampled.\n", "# -beta: Positive integer describing the confidence level required on the bound on the JSR of [2, Prop. 3].\n", "# -beta1: Positive integer describing the confidence level required on the bound of the maximal norm of the modes\n", "# of the system.\n", "# -tol: Positive integer describing the accuracy of the bisection on gamma when solving Problem (1).\n", "# -d: Positive integer describing the degree of the SOS Lyapunov-like function.\n", "# -x0: Float matrix of dimension nxN containing N initial states drawn randomly and uniformly on the sphere.\n", "# -x: Float matrix of dimension nxN containing N final states of randomly drawn trajectories of length l\n", "# associated to the initial states contained in x0.\n", "# -C: Positive integer describing the bound on the norm of P in Problem (1). \n", "#\n", "# Output: -low: lower bound on the JSR of Lemma 3.1 or 3.2\n", "# -upbound: probabilistic bound on the JSR of the system based on sensivity analysis as given in [2, Thm. 3.1] when \n", "# d=1 or [2, Thm. 3.3] otherwise.\n", "# -up: probabilistic bound on the JSR of the system based on chance-constrained optimization as given in [1, Cor.14]\n", "# when d=1, 0 otherwise.\n", "# -upboundbis: probabilistic bound on the JSR of the system based on sensivity analysis as given in [2, Thm. 3.2] when \n", "# d=1, 0 otherwise.\n", "# -Psol: solution of Problem (1).\n", "\n", "function proba_bounds(m,n,l,N,beta,beta1,tol,d,x0,x,C)\n", " U=findmax([norm(x[:,i])/norm(x0[:,i]) for i=1:N])\n", " U=U[1]\n", " D=binomial(n+d-1,d)\n", " #lifting\n", " X=x\n", " X0=x0\n", " x=zeros(D,N)\n", " x0=zeros(D,N)\n", " for i=1:N\n", " x[:,i]=dlift(X[:,i],d,n)\n", " x0[:,i]=dlift(X0[:,i],d,n)\n", " end\n", " ub=U*2\n", " lb=0\n", " gamma=U\n", " #Bisection\n", " iter=1\n", " while ub-lb > tol && iter<100\n", " iter=iter+1\n", " gamma=(ub+lb)/2\n", " solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)\n", " model=Model(solver)\n", " @variable(model, P[1:D, 1:D] in PSDCone())\n", " @SDconstraint(model, P >= Matrix(I,D,D))\n", " @objective(model, Max, 0)\n", " for i=1:N\n", " @constraint(model, x[:,i]'*P*x[:,i] <= gamma^(2*l*d)*x0[:,i]'*P*x0[:,i])\n", " end\n", " @constraint(model, Symmetric(C*Matrix(I,D,D)-P) in PSDCone())\n", " JuMP.optimize!(model)\n", " if termination_status(model) == MOI.OPTIMAL\n", " ub=gamma\n", " else\n", " lb=gamma\n", " end\n", " end\n", " gamma=ub\n", " d1=D*(D+1)/2\n", " #Tie-breaking rule\n", " solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)\n", " model=Model(solver)\n", " @variable(model, P[1:D, 1:D] in PSDCone())\n", " @variable(model, t >=0)\n", " @objective(model, Min, t)\n", " for i=1:N\n", " @constraint(model, x[:,i]'*P*x[:,i] <= gamma^(2*l*d)*x0[:,i]'*P*x0[:,i])\n", " end\n", " @SDconstraint(model, P >= Matrix(I,D,D))\n", " @constraint(model, Symmetric(t*Matrix(I,D,D)-P) in PSDCone())\n", " @constraint(model, Symmetric(C*Matrix(I,D,D)-P) in PSDCone())\n", " optimize!(model)\n", " Psol=value.(P)\n", " #Derivations of the upper bounds\n", " kappa=sqrt(eigmax(Psol)/eigmin(Psol))\n", " kappa2=sqrt(det(Psol)/eigmin(Psol)^D)\n", " epsilon=m^l*(1-((1-beta)/d1)^(1/N))\n", " if epsilon<0.5\n", " delta=sqrt(1-inc_beta_lower_reg_inv((n-1)/2.0,0.5,2*epsilon))\n", " else\n", " delta=0\n", " end\n", " lambda=normmax(m,n,l,N,beta1,X,X0)\n", " if d<2\n", " if N<=d1\n", " epsi=1\n", " else\n", " epsi=1-inc_beta_lower_reg_inv(convert(AbstractFloat, N-d1),convert(AbstractFloat, 4),1-beta+1-beta1)\n", " end\n", " if kappa2*epsi*m^l <1\n", " delt=sqrt(1-inc_beta_lower_reg_inv((n-1)/2.0,0.5,epsi*m^l*kappa2))\n", " elseif abs((1-epsi*m^l)*sqrt(det(Psol)/eigmax(Psol)^n))<0.00001\n", " delt=10^(-10)\n", " elseif abs((1-epsi*m^l)*sqrt(det(Psol)/eigmax(Psol)^n))<0.001\n", " delt=10^(-5)\n", " elseif 1-((1-epsi*m^l)*sqrt(det(Psol)/eigmax(Psol)^n))>1\n", " delt=0\n", " else\n", " delt=sqrt(1-inc_beta_lower_reg_inv((n-1)/2.0,0.5,1-(1-epsi*m^l)*sqrt(det(Psol)/eigmax(Psol)^n)))\n", " end\n", " upold=gamma/delt^(1/l)\n", " upbound=(gamma^l+((gamma^l+lambda)*sqrt(2*(1-delta)))*kappa)^(1/l)\n", " upbounbis=(gamma^(2*l)+(2*(gamma^(2*l)+lambda^2)*sqrt(2*(1-delta))+(-gamma^(2*l)+lambda^2)*2*(1-delta))*kappa^2)^(1/(2*l))\n", " else\n", " upold=0\n", " delt=0\n", " upbounbis=0\n", " upbound=(gamma^(d*l)+((gamma^(d*l)+lambda^d)*sqrt(D)*((1+sqrt(2-2*delta))^d-1)-(gamma^(d*l)+lambda^d)*(sqrt(D)-1)*(sqrt(2-2*delta))^d)*kappa)^(1/(d*l))\n", " end\n", " low=gamma/min(D,m)^(1/(2*d*l))\n", " return (low,upbound, upold,Psol,upbounbis)\n", "end\n", "# function normmax(m,n,l,N,beta1,x,x0)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function computes a probabilistic upper bound on the largest norm of the modes of a system of m modes \n", "# of dimension nxn, based on N trajectories of length l of initial state x0_i and final state x_i for i=1...N.\n", "# The upper bound hold with confidence at least beta1. The computed bound is the one of [2, Proposition 3.2]\n", "#\n", "# Inputs: -m: Positive integer describing the number of modes of the system.\n", "# -n: Positive integer describing the dimension of the modes of the system.\n", "# -l: Positive integer describing the length of the trajectories.\n", "# -N: Positive integer describing the number of trajectories sampled.\n", "# -beta1: Positive integer describing the confidence level required on the bound \n", "# -x0: Float matrix of dimension nxN containing N initial states drawn randomly and uniformly on the sphere.\n", "# -x: Float matrix of dimension nxN containing N final states of randomly drawn trajectories of length l\n", "# associated to the initial states contained in x0.\n", "#\n", "# Output: -lambda*: probabilistic bound largest norm of the modes of the system\n", "\n", "function normmax(m,n,l,N,beta1,x,x0)\n", " M=findmax([norm(x[:,i])/norm(x0[:,i]) for i=1:N])\n", " lambda=M[1]\n", " epsilon1=1-(1-beta1)^(1/N)\n", " delta1=sqrt(inc_beta_upper_reg_inv((n-1)/2.0,0.5,epsilon1*m^l/2))\n", " return (lambda/(delta1))\n", "end" ] }, { "cell_type": "code", "execution_count": 6, "id": "convertible-tunisia", "metadata": {}, "outputs": [ { "ename": "LoadError", "evalue": "\u001b[91mUndefVarError: p not defined\u001b[39m", "output_type": "error", "traceback": [ "\u001b[91mUndefVarError: p not defined\u001b[39m", "", "Stacktrace:", " [1] top-level scope at In[6]:79", " [2] include_string(::Function, ::Module, ::String, ::String) at .\\loading.jl:1091" ] } ], "source": [ "# function comparison(N,Ntrials,m,n,l,d,beta1,beta,tol,A,C)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function analyzes the evolution of probabilistic upper bounds and deterministic lower bounds on the\n", "# JSR of a system A of m modes of dimension nxn, based on N randomly sampled trajectories of length l, \n", "# starting from uniformly, randomly and independently initial states sampled on the sphere, with N. \n", "# The upper bounds hold with confidence at least beta1+beta-1. They are computed using Lyapunov-like \n", "# functions of degree d, for each value of d contained in the vector d. For each value of N contained\n", "# in the vector N and each value of d, the average of the bounds over Ntrials runs is computed. The\n", "# upper bounds are based on [1, Cor. 14], [2, Thm. 3.1], and [2,Thm. 3.2] when d=1; on [2, Thm. 3.3] \n", "# otherwise. The lower bounds are based on [2, Lem. 3.1] when d=1; on [2, Lem. 3.2] otherwise. \n", "# The true value of the JSR of the system is also computed.\n", "#\n", "# For each value of N and d, the function calls the function proba_bounds, which computes an approximate CQLF \n", "# or SOS Lyapunov function of degree 2d by solving the problem \n", "# min gamma s.t. (x_i)^[d]' P (x_i)^[d] ≤ gamma^(2dl) (x0_i)^[d]' P (x0_i)^[d] for i=1...N; P≥I; ||P||≤ C (1)\n", "# where x^[d] is the d-lift of x. The problem is solved by bisection over gamma, with an accuracy tol.\n", "# If d=2, the fonction computes the probabilistic bounds given in [1, Cor. 14] (based on a chance-constrained\n", "# approach), [2, Thm. 3.1] and [2, Thm. 3.2] (based on sensivity analysis). Else, the function computes the upper bound given in \n", "# [2, Thm. 3.3] (based on sensivity analysis).\n", "# For the bounds of [2], two sampled problems have to be solved: (1) and the computation of the maximal norm of the \n", "# modes of the system. beta_1 is the confidence level of the probabilistic upper bound on the maximal norm of the \n", "# modes, while beta is the confidence level on the probabilistic bound on the JSR of the system based on Problem (1)\n", "# and involving the maximal norm of the modes, described in [2, Prop. 3.2]. For all bounds, the total confidence\n", "# level is fixed to beta+beta1-1.\n", "# The default values of the arguments give rise to the first experiment of [2,Chapter 4]\n", "#\n", "# Inputs: -N: Array of positive integers containing all the number of trajectories to be sampled.\n", "# -Ntrials: Positive integer describing the number of trials realised for each value of N and d and over which\n", "# the average of the bounds is computed.\n", "# -m: Positive integer describing the number of modes of the system.\n", "# -n: Positive integer describing the dimension of the modes of the system.\n", "# -l: Positive integer describing the length of the trajectories.\n", "# -d: Array of positive integers containing all the values of degrees of the SOS Lyapunov-like function to be tested. \n", "# -beta: Positive integer describing the confidence level required on the bound on the JSR of [2, Prop. 3].\n", "# -beta1: Positive integer describing the confidence level required on the bound of the maximal norm of the modes\n", "# of the system.\n", "# -tol: Positive integer describing the accuracy of the bisection on gamma when solving Problem (1).\n", "# -A: Matrix containing m float matrices of dimension nxn: the m modes A_1, ..., A_m of dimension nxn of the system.\n", "# -C: Positive integer describing the bound on the norm of P in Problem (1). \n", "#\n", "# Output: -jsr: value of the true JSR of the system\n", "# -up: (size(d))x(size(N)) float matrix. up[1,:] contains the average of the bound of [1,Cor. 14] over Ntrials runs, for each value of N\n", "# -upnew: (size(d))x(size(N)) float matrix. upnew[1,:] contains the average of the bound of [2,Thm. 3.1] over Ntrials runs, for each value of N.\n", "# upnew[i,:] contains the average of the bound of [2,Thm. 3.3] over Ntrials runs, for each value of N, for a degree d[i].\n", "# -low: (size(d))x(size(N)) float matrix. low[i,:] contains the average of the lower bound over Ntrials runs, for each value of N and for a function of degree 2d[i].\n", "\n", "\n", "function comparison(N=append!([300+600*i for i=0:10],[10000]),Ntrials=10,m=2,n=2,l=1,d=[1,2],beta1=0.95,beta=0.95,tol=0.001,A=[[ 1 0;\n", " 1 0 ],\n", "[0 1 ;\n", " 0 -1]],C=20)\n", " p=length(N)\n", " D=length(d)\n", " upnew=zeros(D,p)\n", " up=zeros(D,p)\n", " low=zeros(D,p)\n", " for i=1:Ntrials\n", " for k=1:p\n", " x,x0=generatetrajectories(m,n,l,A,N[k])\n", " for j=1:D\n", " (lo,unew,u,Psol,ubis)=proba_bounds(m,n,l,N[k],beta,beta1,tol,d[j],x0,x,C)\n", " upnew[j,k]+=unew\n", " up[j,k]+=u\n", " low[j,k]+=lo\n", " end\n", " end\n", " end\n", " upnew=upnew/Ntrials\n", " up=up/Ntrials\n", " low=low/Ntrials\n", " jsr=trueJSR(A)\n", " return (jsr,up,upnew,low)\n", "end\n", "#To plot the first experiment of [2, Chapter 4.1]\n", "begin\n", " Plots.plot(N,jsr*ones(1,p)', line = (:dash), label=\"JSR\",linewidth = 2,framestyle=:box)\n", " Plots.plot!(N,up[1,:],label=\"Quadratic case, [BJW21]\",linewidth = 2)\n", " Plots.plot!(N,upnew[1,:],label=\"Quadratic case, Theorem 3.1\",linewidth = 2)\n", " Plots.plot!(N,upnew[2,:],label=\"SOS case, d=2\",linewidth = 2)\n", " lens!([1, 50], [0,25], inset = (1, bbox(0.2, 0.4, 0.7, 0.3)), xticks = 1:5:50, yticks = 0:5:25, subplot = 2, legend=false, framestyle=:box)\n", " ylims!((0.5,5), subplot=1)\n", " xlabel!(\"Number of samples\", subplot=1)\n", " ylabel!(\"Upper bounds on the JSR\", subplot=1) \n", " Plots.savefig(\"JSRBounds.png\")\n", "end\n", "\n", "#To plot the experiment of [2,Chapter 4.5]\n", " begin\n", " Plots.plot(N,jsr*ones(1,p)', line = (:dash), label=\"JSR\",legend=:topleft)\n", " Plots.plot!(N,low[1,:],label=\"Quadratic case\",legend=:topleft)\n", " Plots.plot!(N,low[2,:],label=\"SOS case, d=2\",legend=:topleft)\n", " Plots.plot!(N,low[3,:],label=\"SOS case, d=4\",legend=:topleft)\n", " xlabel!(\"Number of samples\",legend=:topleft)\n", " ylabel!(\"Lower bounds on the JSR\",legend=:topleft)\n", " Plots.savefig(\"JSRLowerBounds.png\")\n", " end" ] }, { "cell_type": "code", "execution_count": 20, "id": "fourth-greece", "metadata": {}, "outputs": [ { "ename": "LoadError", "evalue": "\u001b[91mUndefVarError: p not defined\u001b[39m", "output_type": "error", "traceback": [ "\u001b[91mUndefVarError: p not defined\u001b[39m", "", "Stacktrace:", " [1] top-level scope at .\\In[20]:81", " [2] include_string(::Function, ::Module, ::String, ::String) at .\\loading.jl:1091" ] } ], "source": [ "# function trace(N,Ntrials,m,n,l,d,beta1,beta,tol,A,C)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function analyzes the evolution of probabilistic upper bounds and deterministic lower bounds on the\n", "# JSR of a system A of m modes of dimension nxn, based on N randomly sampled trajectories of length l for various values of l, \n", "# starting from uniformly, randomly and independently initial states sampled on the sphere, with N. \n", "# The upper bounds hold with confidence at least beta1+beta-1. They are computed using Lyapunov-like \n", "# functions of degree d. For each value of N contained\n", "# in the vector N and each value of l, the average of the bounds over Ntrials runs is computed. The\n", "# upper bounds are based on [1, Cor. 14], [2, Thm. 3.1], and [2,Thm. 3.2] when d=1; on [2, Thm. 3.3] \n", "# otherwise. The lower bounds are based on [2, Lem. 3.1] when d=1; on [2, Lem. 3.2] otherwise. \n", "# The true value of the JSR of the system is also computed.\n", "#\n", "# For each value of N and l, the function calls the function proba_bounds, which computes an approximate CQLF \n", "# or SOS Lyapunov function of degree 2d by solving the problem \n", "# min gamma s.t. (x_i)^[d]' P (x_i)^[d] ≤ gamma^(2dl) (x0_i)^[d]' P (x0_i)^[d] for i=1...N; P≥I; ||P||≤ C (1)\n", "# where x^[d] is the d-lift of x. The problem is solved by bisection over gamma, with an accuracy tol.\n", "# If d=2, the fonction computes the probabilistic bounds given in [1, Cor. 14] (based on a chance-constrained\n", "# approach), [2, Thm. 3.1] and [2, Thm. 3.2] (based on sensivity analysis). Else, the function computes the upper bound given in \n", "# [2, Thm. 3.3] (based on sensivity analysis).\n", "# For the bounds of [2], two sampled problems have to be solved: (1) and the computation of the maximal norm of the \n", "# modes of the system. beta_1 is the confidence level of the probabilistic upper bound on the maximal norm of the \n", "# modes, while beta is the confidence level on the probabilistic bound on the JSR of the system based on Problem (1)\n", "# and involving the maximal norm of the modes, described in [2, Prop. 3.2]. For all bounds, the total confidence\n", "# level is fixed to beta+beta1-1.\n", "# The default values of the arguments give rise to the first experiment of [2,Chapter 4]\n", "#\n", "# Inputs: -N: Array of positive integers containing all the number of trajectories to be sampled.\n", "# -Ntrials: Positive integer describing the number of trials realised for each value of N and d and over which\n", "# the average of the bounds is computed.\n", "# -m: Positive integer describing the number of modes of the system.\n", "# -n: Positive integer describing the dimension of the modes of the system.\n", "# -l: Array of positiv integers describing the length of the trajectories.\n", "# -d: Positive integer describing the degree of the SOS Lyapunov-like function to be tested. \n", "# -beta: Positive integer describing the confidence level required on the bound on the JSR of [2, Prop. 3].\n", "# -beta1: Positive integer describing the confidence level required on the bound of the maximal norm of the modes\n", "# of the system.\n", "# -tol: Positive integer describing the accuracy of the bisection on gamma when solving Problem (1).\n", "# -A: Matrix containing m float matrices of dimension nxn: the m modes A_1, ..., A_m of dimension nxn of the system.\n", "# -C: Positive integer describing the bound on the norm of P in Problem (1). \n", "#\n", "# Output: -up: (size(l))x(size(N)) float matrix. up[i,:] contains the average of the bound of [1,Cor. 14] over Ntrials runs, for each value of N, for l[i]\n", "# -upnew: (size(l))x(size(N)) float matrix. upnew[i,:] contains the average of the bound of [2,Thm. 3.1] over Ntrials runs, for each value of N, for l[i].\n", "\n", "function trace(Ntrials=10,N=append!([200+300*i for i=0:10],[4000]),d=1,l=[1,3,5],A=generatesystem(3,3),m=3,n=3,beta1=0.95,beta=0.95,tol=0.001,C=100)\n", " p=length(N)\n", " D=length(l)\n", " upnew=zeros(D,p)\n", " up=zeros(D,p)\n", " for i=1:Ntrials\n", " for j=1:D\n", " x0=zeros(n,N[p])\n", " x=zeros(n,N[p])\n", " for k=1:p\n", " if k==1\n", " X0,X=generatetrajectories(m,n,l[j],A,N[k])\n", " x0[:,1:N[k]]=X0\n", " x[:,1:N[k]]=X\n", " else\n", " X0,X=generatetrajectories(m,n,l[j],A,N[k]-N[k-1])\n", " x0[:,1+N[k-1]:N[k]]=X0\n", " x[:,1+N[k-1]:N[k]]=X\n", " end\n", " (lo,unew, u, Psol, upboundbis)=proba_bounds(m,n,l[j],N[k],beta,beta1,tol,d,x0,x,C)\n", " upnew[j,k]+=unew\n", " up[j,k]+=u\n", " end\n", " end\n", " end\n", " upnew=upnew/Ntrials\n", " up=up/Ntrials\n", " jsr=trueJSR(A)\n", " \n", " return(up,upnew)\n", "end\n", "#To plot the experiment of [2, Chapter 4.2]\n", " begin\n", " plot(N,jsr*ones(1,p)', line = (:dash), label=\"JSR\",linewidth = 2, framestyle=:box)\n", " for j=1:3\n", " plot!(N,upnew[j,:], label=\"l=\" * string(l[j]),linewidth = 2)\n", " end\n", " xlabel!(\"Number of samples\")\n", " ylabel!(\"Upper bounds on the JSR\")\n", " ylims!((0,5))\n", " savefig(\"JSRtraced2.png\")\n", " end" ] }, { "cell_type": "code", "execution_count": 8, "id": "uniform-symbol", "metadata": {}, "outputs": [ { "ename": "LoadError", "evalue": "\u001b[91mUndefVarError: up not defined\u001b[39m", "output_type": "error", "traceback": [ "\u001b[91mUndefVarError: up not defined\u001b[39m", "", "Stacktrace:", " [1] top-level scope at .\\In[8]:84", " [2] include_string(::Function, ::Module, ::String, ::String) at .\\loading.jl:1091" ] } ], "source": [ "# function convergence(Ntrials,N,d,m,n,l,beta1,beta,tol,A,C)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function analyzes the convergence of probabilistic upper bounds on the\n", "# JSR of systems A of various dimensions, based on N randomly sampled trajectories of length l, \n", "# starting from uniformly, randomly and independently initial states sampled on the sphere, when N increases. \n", "# The upper bounds hold with confidence at least beta1+beta-1. They are computed using Lyapunov-like \n", "# functions of degree d. For each combination of m and n, a system is generated using the function generatesystem.\n", "# For each value of N contained in the vector N and each value of m,n;\n", "# the average of the bounds over Ntrials runs is computed. The\n", "# upper bounds are based on [1, Cor. 14], [2, Thm. 3.1], and [2,Thm. 3.2] when d=1; on [2, Thm. 3.3] \n", "# otherwise. The ratio of the bounds on the true value of the JSR of the system is returned.\n", "#\n", "# For each value of N and d, the function calls the function proba_bounds, which computes an approximate CQLF \n", "# or SOS Lyapunov function of degree 2d by solving the problem \n", "# min gamma s.t. (x_i)^[d]' P (x_i)^[d] ≤ gamma^(2dl) (x0_i)^[d]' P (x0_i)^[d] for i=1...N; P≥I; ||P||≤ C (1)\n", "# where x^[d] is the d-lift of x. The problem is solved by bisection over gamma, with an accuracy tol.\n", "# If d=2, the fonction computes the probabilistic bounds given in [1, Cor. 14] (based on a chance-constrained\n", "# approach), [2, Thm. 3.1] and [2, Thm. 3.2] (based on sensivity analysis). Else, the function computes the upper bound given in \n", "# [2, Thm. 3.3] (based on sensivity analysis).\n", "# For the bounds of [2], two sampled problems have to be solved: (1) and the computation of the maximal norm of the \n", "# modes of the system. beta_1 is the confidence level of the probabilistic upper bound on the maximal norm of the \n", "# modes, while beta is the confidence level on the probabilistic bound on the JSR of the system based on Problem (1)\n", "# and involving the maximal norm of the modes, described in [2, Prop. 3.2]. For all bounds, the total confidence\n", "# level is fixed to beta+beta1-1.\n", "# The default values of the arguments give rise to the first experiment of [2,Chapter 4.3]\n", "#\n", "# Inputs: -N: Array of positive integers containing all the number of trajectories to be sampled.\n", "# -Ntrials: Positive integer describing the number of trials realised for each value of N and d and over which\n", "# the average of the bounds is computed.\n", "# -m: Array of positive integers describing the number of modes of the system.\n", "# -n: Array of positive integers describing the dimension of the modes of the system.\n", "# -l: Positive integer describing the length of the trajectories.\n", "# -d: Positive integers describing the degree of the SOS Lyapunov-like function to be tested. \n", "# -beta: Positive integer describing the confidence level required on the bound on the JSR of [2, Prop. 3].\n", "# -beta1: Positive integer describing the confidence level required on the bound of the maximal norm of the modes\n", "# of the system.\n", "# -tol: Positive integer describing the accuracy of the bisection on gamma when solving Problem (1).\n", "# -A: Matrix containing m float matrices of dimension nxn: the m modes A_1, ..., A_m of dimension nxn of the system.\n", "# -C: Positive integer describing the bound on the norm of P in Problem (1). \n", "#\n", "# Output: -up: (size(n))x(size(N)) float matrix. up[i,:] contains the average of the ratio of the bound of [1,Cor. 14] over the true value of the JSR over Ntrials runs, for each value of N, for a system of dimension m[i],n[i].\n", "# -upnew: (size(d))x(size(N)) float matrix. upnew[i,:] contains the average of the ratio of the bound of [2,Thm. 3.1] over the true value of the JSR over Ntrials runs, for each value of N, for a system of dimension m[i],n[i].\n", "# -upnewbis: (size(d))x(size(N)) float matrix. upnewbis[i,:] contains the average of the ratio of the bound of [2,Thm. 3.2] over the true value of the JSR over Ntrials runs, for each value of N, for a system of dimension m[i],n[i].\n", "\n", "function convergence(Ntrials=10,N=append!([100+300*i for i=0:10],[4000]),d=1,m=[5,5,7],n=[3,5,3],l=1,beta1=0.95,beta=0.95,tol=0.001,C=100)\n", " p=length(N)\n", " D=length(m)\n", " upnew=zeros(D,p)\n", " upnewbis=zeros(D,p)\n", " up=zeros(D,p)\n", " jsr=zeros(D,1)\n", " for i=1:Ntrials\n", " for j=1:D\n", " A=generatesystem(m[j],n[j])\n", " jsr[j]=trueJSR(A)\n", " x0=zeros(n[j],N[p])\n", " x=zeros(n[j],N[p])\n", " for k=1:p\n", " if k==1\n", " X0,X=generatetrajectories(m[j],n[j],l,A,N[k])\n", " x0[:,1:N[k]]=X0\n", " x[:,1:N[k]]=X\n", " else\n", " X0,X=generatetrajectories(m[j],n[j],l,A,N[k]-N[k-1])\n", " x0[:,1+N[k-1]:N[k]]=X0\n", " x[:,1+N[k-1]:N[k]]=X\n", " end\n", " (lo,unew,u,Psol,upboundbis)=proba_bounds(m[j],n[j],l,N[k],beta,beta1,tol,d,x0,x,C)\n", " upnew[j,k]+=unew/jsr[j]\n", " upnewbis[j,k]+=upboundbis/jsr[j]\n", " up[j,k]+=u/jsr[j]\n", " end\n", " end\n", " end\n", " return(up/Ntrials,upnew/Ntrials,upnewbis/Ntrials)\n", "end\n", "#To plot the experiment of [2, Chapter 4.3]\n", " begin\n", " Plots.plot(N,up[1,:], label=\"[BJW21], n=\" * string(n[1]) * \", m=\" * string(m[1]),linewidth = 2, framestyle=:box)\n", " Plots.plot!(N,upnew[1,:], label=\"Theorem 3.1, n=\" * string(n[1]) * \", m=\" * string(m[1]),linewidth = 2)\n", " for j=2:2\n", " Plots.plot!(N,up[j,:], label=\"[BJW21], n=\" * string(n[j]) * \", m=\" * string(m[j]),linewidth = 2)\n", " Plots.plot!(N,upnew[j,:], label=\"Theorem 3.1, n=\" * string(n[j]) * \", m=\" * string(m[j]),linewidth = 2)\n", " end\n", " xlabel!(\"Number of samples\")\n", " ylabel!(\"Black-box bound/White-box bound\")\n", " ylims!((0,10))\n", " Plots.savefig(\"JSRconvergence.png\")\n", " end\n", " begin\n", " Plots.plot(N,up[1,:], label=\"[BJW21], n=\" * string(n[1]) * \", m=\" * string(m[1]),linewidth = 2, framestyle=:box)\n", " Plots.plot!(N,upnew[1,:], label=\"Theorem 3.1, n=\" * string(n[1]) * \", m=\" * string(m[1]),linewidth = 2)\n", " for j=3:3\n", " Plots.plot!(N,up[j,:], label=\"[BJW21], n=\" * string(n[j]) * \", m=\" * string(m[j]),linewidth = 2)\n", " Plots.plot!(N,upnew[j,:], label=\"Theorem 3.1, n=\" * string(n[j]) * \", m=\" * string(m[j]),linewidth = 2)\n", " end\n", " xlabel!(\"Number of samples\")\n", " ylabel!(\"Black-box bound/White-box bound\")\n", " ylims!((0,10))\n", " Plots.savefig(\"JSRconvergence1.png\")\n", " end\n", "#To plot the experiment of [2, Chapter 4.4]\n", " begin\n", " plot(N,upnew[1,:], label=\"Bound of Theorem 3.1, n=\" * string(n[1]) * \", m=\" * string(m[1]),linewidth = 2, framestyle=:box)\n", " plot!(N,upnewbis[1,:], label=\"Bound of Theorem 3.2, n=\" * string(n[1]) * \", m=\" * string(m[1]),linewidth = 2)\n", " for j=2:2\n", " plot!(N,upnew[j,:], label=\"Bound of Theorem 3.1, n=\" * string(n[j]) * \", m=\" * string(m[j]),linewidth = 2)\n", " plot!(N,upnewbis[j,:], label=\"Bound of Theorem 3.2, n=\" * string(n[j]) * \", m=\" * string(m[j]),linewidth = 2)\n", " end\n", " xlabel!(\"Number of samples\")\n", " ylabel!(\"Upper bound (black-box)/Upper bound (white-box)\")\n", " lens!([0, 400], [1,5], inset = (1, bbox(0.2, 0.4, 0.7, 0.3)))\n", " savefig(\"JSRconvergence.png\")\n", " end" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.5.3", "language": "julia", "name": "julia-1.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.5.3" } }, "nbformat": 4, "nbformat_minor": 5 }