{ "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": 3, "id": "speaking-review", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "truesol (generic function with 1 method)" ] }, "execution_count": 3, "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 \n", "# function truesol(A,m,n,d,C,tol)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function computes an upper bound on the JSR of a system of m modes of dimension nxn, in a white-box way.\n", "# The function computes an approximate CQLF or SOS Lyapunov function of degree 2d by solving the problem \n", "# min gamma s.t. (A_i)^[d]' P (A)^[d] ≤ gamma^(2dl) P for i=1...m; P≥I; ||P||≤ C (1)\n", "# where A^[d] is the d-lift of A. The problem is solved by bisection over gamma, with an accuracy tol.\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", "# -A: 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", "# -C: Positive integer describing the bound on the norm of P in Problem (1). \n", "#\n", "#Output: -Psol, solution to Problem (1)\n", "function truesol(A,m,n,d,C,tol)\n", " l=1\n", " A1=[zeros(n,n) for i=1:m]\n", " ub=0\n", " lb=0\n", " for i=1:m\n", " A1[i]=dliftmat(A[i],d,n)\n", " if opnorm(A1[i],2)>ub\n", " ub=opnorm(A1[i],2)\n", " end\n", " end\n", " gamma=0\n", " D=binomial(n+d-1,d)\n", " ub=ub*2\n", " while ub-lb > tol\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:m\n", " @SDconstraint(model, A1[i]'*P*A1[i] <= gamma^(2*l*d)*P)\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", " 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:m\n", " @SDconstraint(model, A1[i]'*P*A1[i] <= gamma^(2*l*d)*P)\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", " return (Psol)\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": 5, "id": "stuffed-northern", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "inversedlift (generic function with 1 method)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# function inversedlift(Ld,d,n)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function computes the least square solution L, of dimension nxn,\n", "#to the following problem: min ||L^d-Ld||, where L^d is the d-lift of L and Ld is a\n", "#DxD matrix.\n", "#\n", "# Inputs: -Ld: Float matrix of dimension D=binomial(n+d-1,d).\n", "# -d: Positive integer describing the dimension of the d-lift.\n", "# -n: Positive integer describing the dimension of the matrix L to compute. \n", "#\n", "# Output: -L: Float matrix of dimension n, least square solution to min ||L^d-Ld||.\n", "\n", "function inversedlift(Ld,d,n)\n", " D=binomial(n+d-1,d)\n", " A=zeros(1,n)\n", " B=zeros(1,n)\n", " C=zeros(1,n)\n", " E=zeros(1,n)\n", " M=zeros(d,d)\n", " alpha=weakcompo(d,n)\n", " model = Model(Ipopt.Optimizer)\n", " @variable(model,aux[1:D,1:D])\n", " @variable(model, L[1:n,1:n])\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", " my_expr=@NLexpression(model,sum(prod(L[convert(Int64,C[k]),convert(Int64,E[y[p][k]])] for k=1:d) for p=1:length(y)))\n", " @NLconstraint(model,aux[i,j]==my_expr/den)\n", " set_start_value(aux[i,j], Ld[i,j])\n", " end\n", " end\n", " for i=1:n\n", " for j=1:n\n", " set_start_value(L[i,j],abs(Ld[1+D-convert(Int64,binomial(n+d-i,d)),1+D-convert(Int64,binomial(n+d-j,d))])^(1/d))\n", " end\n", " end \n", " @NLobjective(model,Min,sum(sum((aux[i,j]-Ld[i,j])^2 for i=1:D) for j=1:D))\n", " JuMP.optimize!(model)\n", " return(value.(L))\n", "end" ] }, { "cell_type": "code", "execution_count": 6, "id": "gorgeous-animation", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "normmax (generic function with 1 method)" ] }, "execution_count": 6, "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", "# 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: -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", "# -gamma, Psol: solution to Problem (1).\n", "# -X0,X=x0,x\n", "# -lambda: probabilistic bound on the largest norm of the modes of the system \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", " 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", " 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", " 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", " up=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", " up=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", " return (upbound, up, gamma,X,X0,Psol,lambda,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 on the largest norm of the modes of the system\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": 7, "id": "rational-brook", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "normmaxbis (generic function with 1 method)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# function proba_bounds_offline_sampling(m,n,l,N,beta,beta1,tol,d,X0,X,C,xx0,xx,N1,L)\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 L A Linv 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, and on N1 trajectories of \n", "# length l of initial state L*xx0_i and final state L*xx_i (non uniform sampling), that is by solving Opt_k(w_1+...+w_k).\n", "# The upper bounds hold 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+N1; 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=1, the fonction computes the probabilistic bounds given in [2, Thm. 5.1] and [2,Thm. 5.2]. Else, the function computes the upper bound given in \n", "# [2, Thm.6.1].\n", "# 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 uniformly in the system.\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: X_i=LA_jLinv X0_i.\n", "# -C: Positive integer describing the bound on the norm of P in Problem (1). \n", "# -xx0: Float matrix of dimension nxN1 containing N1 initial states drawn randomly and uniformly on the sphere.\n", "# -xx: Float matrix of dimension nxN1 containing N1 final states of randomly drawn trajectories of length l\n", "# associated to the initial states contained in xx0: xx_i=A_j xx0_i.\n", "# -N1: Positive integer describing the number of trajectories sampled non uniformly in the system.\n", "# -L: matrix determining the coordinate system in which X0,X is uniformly sampled. \n", "#\n", "# Output: -upbound: probabilistic bound on the JSR of the system based on sensivity analysis as given in [2, Thm. 5.1] when \n", "# d=1 or [2, Thm. 6.1] 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", "# -gamma,Psol: solution of Problem (1).\n", "# -X0=union of X0 and Lxx0\n", "# -X=union of X and Lxx\n", "# -lambda: probabilistic bound on the JSR of the system based on sensivity analysis as given in [2, Thm. 5.2] when \n", "# d=1, 0 otherwise.\n", "\n", "function proba_bounds_offline_sampling(m,n,l,N,beta,beta1,tol,d,X0,X,C,xx0,xx,N1,L)\n", " D=binomial(n+d-1,d)\n", " x0=zeros(n,N+N1)\n", " x=zeros(n,N+N1)\n", " x0[:,1:N]=X0\n", " x[:,1:N]=X\n", " for i=1:N1\n", " x0[:,i+N]=L*xx0[:,i]\n", " x[:,i+N]=L*xx[:,i]\n", " end\n", " U=findmax([norm(x[:,i])/norm(x0[:,i]) for i=1:N+N1])\n", " U=U[1]\n", " #lifting\n", " X=x\n", " X0=x0\n", " x=zeros(D,N+N1)\n", " x0=zeros(D,N+N1)\n", " for i=1:N+N1\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", " 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+N1\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", " 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+N1\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", " 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=normmaxbis(m,n,l,N,N1,beta1,X,X0)\n", " if d<2\n", " if N-d1>0\n", " epsi=1-inc_beta_lower_reg_inv(convert(AbstractFloat, N-d1),convert(AbstractFloat, 4),1-beta+1-beta1)\n", " else\n", " epsi=1\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", " up=gamma/delt^(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", " upbound=(gamma^l+((gamma^l+lambda)*sqrt(2*(1-delta)))*kappa)^(1/l)\n", " else\n", " up=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", " return (upbound, gamma,Psol, lambda,X,X0,upbounbis)\n", "end\n", "function normmaxbis(m,n,l,N, N1,beta1,x,x0)\n", " M=findmax([norm(x[:,i])/norm(x0[:,i]) for i=1:N+N1])\n", " lambda=M[1]\n", " epsilon1=1-(1-beta1)^(1/N)\n", " if epsilon1*m^l/2>0.99\n", " delta1=10^(-10)\n", " else\n", " delta1=sqrt(1-inc_beta_lower_reg_inv((n-1)/2.0,0.5,epsilon1*m^l/2))\n", " end\n", " return (lambda/(delta1))\n", "end" ] }, { "cell_type": "code", "execution_count": 8, "id": "altered-specialist", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "procedure_online_sampling (generic function with 1 method)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# function procedure_online_sampling(m,n,l,N,N1,Ntot,alpha,tol,d,X0,X,C,xx0,xx,L,Pfinal)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function computes the number of samples required to reach convergence as defined in [2, Procedure Online]. The\n", "# maximum number of trajectories of length l is of Ntot. We start the iterations with N1 trajectories of initial state L*xx0_i\n", "# and final state L*xx_i (non uniform sampling). Then, we add uniformly sampled trajectories (X0,X).\n", "# For i=0,1,2,..., 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...N1+N*i; 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", "# We stop the iterations when the norm of the difference between two successive solutions P is smaller than alpha.\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 increment in trajectories at each iteration.\n", "# -Ntot: Positive integer describing the maximal number of trajectories sampled uniformly in the system.\n", "# -alpha: accuracy on the convergence of the solution P to Problem (1)\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 nxNtot-N1 containing Ntot-N1 initial states drawn randomly and uniformly on the sphere.\n", "# -X: Float matrix of dimension nxNtot-N1 containing Ntot-N1 final states of randomly drawn trajectories of length l\n", "# associated to the initial states contained in X0: X_i=LA_jLinv X0_i.\n", "# -C: Positive integer describing the bound on the norm of P in Problem (1). \n", "# -Pfinal: solution to truesol(A)\n", "# -xx0: Float matrix of dimension nxN1 containing N1 initial states drawn randomly and uniformly on the sphere.\n", "# -xx: Float matrix of dimension nxN1 containing N1 final states of randomly drawn trajectories of length l\n", "# associated to the initial states contained in xx0: xx_i=A_j xx0_i.\n", "# -N1: Positive integer describing the number of trajectories sampled non uniformly in the system.\n", "# -L: matrix determining the coordinate system in which X0,X is uniformly sampled. \n", "#\n", "# Output: -gamma,Psol: solution to Problem (1).\n", "# -X0=union of xx0 and the used samples of X0\n", "# -X=union of xx and the used samples of X.\n", "# -Nused: number of trajectories used to reach convergence on P\n", "# -condition: Evolution of the convergence of P, and of its convergence to Pfinal\n", "\n", "function procedure_online_sampling(m,n,l,N,N1,Ntot,alpha,tol,d,X0,X,C,xx0,xx,L,Pfinal)\n", " condition=zeros(2,convert(Int64, ceil((Ntot-N1)/N))+100)\n", " U=0\n", " D=binomial(d+n-1,d)\n", " gamma=0\n", " upbound=0\n", " up=0\n", " lambda=0\n", " x0=zeros(n,Ntot)\n", " x=zeros(n,Ntot)\n", " for i=1:N1\n", " x0[:,i]=L*xx0[:,i]\n", " x[:,i]=L*xx[:,i]\n", " end\n", " Psol=zeros(n,n)\n", " adapt=0\n", " nono=100\n", " x0[:,1+N*adapt+N1:Ntot]=X0[:,1:Ntot-N1] \n", " x[:,1+N*adapt+N1:Ntot]=X[:,1:Ntot-N1]\n", " XX0=x0\n", " XX=x\n", " x=zeros(D,Ntot)\n", " x0=zeros(D,Ntot)\n", " for i=1:Ntot\n", " x[:,i]=dlift(XX[:,i],d,n)\n", " x0[:,i]=dlift(XX0[:,i],d,n)\n", " end\n", " N0=0\n", " if N1==0\n", " N0=convert(Int64,n*(n+1)/2+1)\n", " end \n", " U=findmax([norm(x[:,i])/norm(x0[:,i]) for i=1:N1+N0])\n", " U=U[1]\n", " ub=U*2\n", " lb=0\n", " gamma=U\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*adapt+N1+N0\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=n*(n+1)/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", " @variable(model, t >=0)\n", " @objective(model, Min, t)\n", " for i=1:N*adapt+N1+N0\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", " condition[1,1]=opnorm(Psol-Pfinal,2)\n", " kappa=sqrt(eigmax(Psol)/eigmin(Psol))\n", " kappa2=sqrt(det(Psol)/eigmin(Psol)^D)\n", " while nono>alpha && N*(adapt+1) 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*adapt+N1+N0\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=n*(n+1)/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", " @variable(model, t >=0)\n", " @objective(model, Min, t)\n", " for i=1:N*adapt+N1+N0\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", " Psolold=Psol\n", " Psol=value.(P)\n", " nono=opnorm(Psol-Psolold,2)\n", " kappa=sqrt(eigmax(Psol)/eigmin(Psol))\n", " kappa2=sqrt(det(Psol)/eigmin(Psol)^D)\n", " condition[1,adapt+1]=opnorm(Psol-Pfinal,2)\n", " condition[2,adapt+1]=nono\n", " end\n", " Nused=adapt*N+N0\n", " return (gamma,Psol,XX,XX0,Nused, condition[:,1:adapt+1])\n", "end" ] }, { "cell_type": "code", "execution_count": 9, "id": "harmful-second", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "f! (generic function with 1 method)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# function algo_fast(m,n,l,N1,N2,beta,beta1,tol,d,x0,x, C)\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function implements the fast algorithm with N1 samples used in the first phase and N2 samples used in the second phase. computes probabilistic upper bounds on the JSR of a system L A Linv of m modes of dimension nxn, based\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", "# -N1: Positive integer describing the number of trajectories used in the first phase.\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 nxN1+N2 containing initial states drawn randomly and uniformly on the sphere.\n", "# -x: Float matrix of dimension nxN1+N2 containing final states of randomly drawn trajectories of length l\n", "# associated to the initial states contained in X0: X_i=A_j X0_i.\n", "# -C: Positive integer describing the bound on the norm of P in Problem (1). \n", "# -N2: Positive integer describing the number of trajectories used in the second phase.\n", "#\n", "# Output: -up: probabilistic bound on the JSR of the system based on N1+N2 as given in [2, Thm. 3.1]\n", "# -up1: probabilistic bound on the JSR of the system based on the fast algorithm as given in [2, Prop. 5.1] when \n", "function algo_fast(m,n,l,N1,N2,beta,beta1,tol,d,x0,x, C)\n", " X0=x0[:,1:N1]\n", " X=x[:,1:N1]\n", " (u, up, gam,xx,xx0,PP,ubis)=proba_bounds(m,n,l,N1,beta,beta1,tol,d,X0,X,C)\n", " X0=x0[:,1+N1:N1+N2]\n", " X=x[:,1+N1:N1+N2]\n", " (low1, up1, gam1)=proba_bounds_fast(m,n,l,N2,beta,beta1,d,X0,X,PP,gam,N1)\n", " return (up,up1)\n", "end\n", "function proba_bounds_fast(m,n,l,N,beta,beta1,d,x0,x,P,gam,N1)\n", " U=findmax([norm(x[:,i])/norm(x0[:,i]) for i=1:N])\n", " U=U[1]\n", " D=n\n", " d1=D*(D+1)/2\n", " L=[(convert(AbstractFloat, x[:,i]'*P*x[:,i])/convert(AbstractFloat, x0[:,i]'*P*x0[:,i]))^(1/(2*l*d)) for i=1:N] \n", " gamma=findmax(L)\n", " gamma=findmax([gamma[1],gam^(1/l)])\n", " gamma=gamma[1]\n", " kappa=sqrt(eigmax(P)/eigmin(P))\n", " kappa2=sqrt(det(P)/eigmin(P)^D)\n", " if N1-d1<11\n", " epsi=1\n", " else\n", " epsi=nlsolve((F,x) ->f!(F,x,beta,beta1,N,N1,d1), [0.001])\n", " epsi=(epsi.zero)[1]\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(P)/eigmax(P)^n))<0.00001\n", " delt=10^(-10)\n", " elseif abs((1-epsi*m^l)*sqrt(det(P)/eigmax(P)^n))<0.001\n", " delt=10^(-5)\n", " elseif 1-((1-epsi*m^l)*sqrt(det(P)/eigmax(P)^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(P)/eigmax(P)^n)))\n", " end\n", " up=gamma/delt^(1/l)\n", " return (gamma/min(D,m)^(1/(2*d*l)),up, gamma)\n", "end\n", "function f!(F, x,beta,beta1,N,N1,d1)\n", " F[1] = 2-beta-beta1-(1-x[1])^N*inc_beta_lower_reg(convert(AbstractFloat, N1-d1),convert(AbstractFloat, 3),1-x[1])\n", "end" ] }, { "cell_type": "code", "execution_count": 10, "id": "destroyed-auction", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "algo_offline_sampling (generic function with 1 method)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# function algo_offline_sampling(m,n,l,N1,N2,beta,beta1,tol,d,x0,x,A, C, phases)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function implements [2, Algorithm 3]. It computes probabilistic bounds on the JSR of a system A of m modes of dimension n. \n", "# It is based on phases+1 phases. In the first phase, one uses N1 samples. In the next phases, one uses N2/phases samples, from the set x0,x.\n", "# For i=0,1,2,..., 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...N1+N2/phases*i; 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", "# 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", "# -N1: Positive integer describing the number of trajectories to use in the first phase.\n", "# -N2: Positive integer describing the number of trajectories to use in the next phases.\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 nxN1+N2 containing initial states drawn randomly and uniformly on the sphere.\n", "# -x: Float matrix of dimension nxN1+N2 containing final states of randomly drawn trajectories of length l\n", "# associated to the initial states contained in X0: X_i=A_j X0_i.\n", "# -C: Positive integer describing the bound on the norm of P in Problem (1). \n", "# -A: system.\n", "# -phases: number of phases (-1) of [2, Algorithm 3]\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", "#\n", "# Output: -u: probabilistic bound on the JSR of the system based on N1+N2 samples as given in [2, Thm. 3.1] when \n", "# d=1 or [2, Thm. 3.3] otherwise.\n", "# -u1: Vector of dimension phases, containing the evolution of the probabilistic bound on the JSR of the\n", "# system with each phase, as given in [2, Thm. 5.1] when d=1 or [2, Thm. 6.1] otherwise.\n", "# -U1: Vector of dimension phases, containing the evolution of the probabilistic bound on the JSR of the\n", "# system with each phase, as given in [2, Thm. 5.2] when d=1 or 0 otherwise.\n", "# -kappa: condition number of the solution Psol to Problem (1) based on N1+N2 samples.\n", "# -kappa1: Vector of dimension phases, containing the evolution of the condition number of the solution Psol to Problem (1)\n", "# with each phase.\n", "\n", "\n", "function algo_offline_sampling(m,n,l,N1,N2,beta,beta1,tol,d,x0,x,A, C, phases)\n", " X0=x0[:,1:N1]\n", " X=x[:,1:N1]\n", " (u, uj, gam,xx,xx0,PP,lam,ubis)=proba_bounds(m,n,l,N1,beta,beta1,tol,d,X0,X,C)\n", " L=cholesky(PP)\n", " L=L.U\n", " Linv=inv(L)\n", " A1=[zeros(n,n) for i=1:m]\n", " for i=1:m\n", " A1[i]=A[i]\n", " end\n", " PP1=[zeros(n,n) for i=1:phases]\n", " u1=zeros(1,phases)\n", " uj1=zeros(1,phases)\n", " gam1=zeros(1,phases)\n", " U1=zeros(1,phases)\n", " for i=1:phases\n", " X=zeros(n,N2)\n", " X0=x0[:,1+N1+(i-1)*N2:N1+i*N2]\n", " for j=1:N2\n", " b=rand((1:m),l)\n", " X[:,j]=Linv*X0[:,j] \n", " for k=1:l\n", " X[:,j]=A1[b[k]]*X[:,j]\n", " end\n", " X[:,j]=L*X[:,j] \n", " end\n", " for j=1:m\n", " A1[j]=L*A1[j]*Linv\n", " end\n", " (u1[i], gam1[i], PP1[i],U1[i],xx,xx0,ubis)=proba_bounds_offline_sampling(m,n,l,N2,beta,beta1,tol,d,X0,X,C,xx0,xx,N1+(i-1)*N2,L)\n", " L=cholesky(PP1[i])\n", " L=L.U\n", " Linv=inv(L)\n", " end\n", " return (u,u1,U1,sqrt(eigmax(PP)/eigmin(PP)),[sqrt(eigmax(PP1[i])/eigmin(PP1[i])^n) for i=1:phases])\n", "end" ] }, { "cell_type": "code", "execution_count": 11, "id": "radical-macintosh", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "algo_online_sampling (generic function with 1 method)" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# function algo_online_sampling(m,n,l,N,N1,beta,beta1,tol,d,x0,x, C, A, phases, alpha)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function implements [2, Algorithm 4]. It computes probabilistic bounds on the JSR of a system A of m modes of dimension n. \n", "# The number of phases used is the minimum between phases and the number decided by [2,Algorithm 4]. \n", "# The total number of samples used is of N1. In each phase, the number of samples used is determined by procedure_online_sampling.\n", "# The increment is fixed to N. The samples are taken from the set x0,X\n", "# For i=0,1,2,..., 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...N1+N2/phases*i; 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", "# 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 increment in trajectories at each iteration of Procedure online.\n", "# -N1: Positive integer describing the total number of trajectories.\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: X_i=A_j X0_i.\n", "# -C: Positive integer describing the bound on the norm of P in Problem (1). \n", "# -A: system\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", "# -alpha: accuracy on the convergence of the solution P to Problem (1)\n", "#\n", "# Output: -u: probabilistic bound on the JSR of the system based on N1 samples as given in [2, Thm. 3.1] when \n", "# d=1 or [2, Thm. 3.3] otherwise.\n", "# -u1: Vector of dimension Nphases, containing the evolution of the probabilistic bound on the JSR of the\n", "# system with each phase, as given in [2, Thm. 5.1] when d=1 or [2, Thm. 6.1] otherwise.\n", "# -U1: Vector of dimension Nphases, containing the evolution of the probabilistic bound on the JSR of the\n", "# system with each phase, as given in [2, Thm. 5.2] when d=1 or 0 otherwise.\n", "# -kappa: condition number of the solution Psol to Problem (1) based on N1 samples.\n", "# -kappa1: Vector of dimension Nphases, containing the evolution of the condition number of the solution Psol to Problem (1)\n", "# with each phase.\n", "# -Nres: Vector of dimension Nphases containing the number of samples used in each phase.\n", "# -Nphases: Number of phases used.\n", "\n", "function algo_online_sampling(m,n,l,N, N1,beta,beta1,tol,d,x0,x, C, A, phases, alpha)\n", " D=binomial(d+n-1,d)\n", " PP1=[zeros(D,D) for i=1:phases+1]\n", " Nphases=1\n", " X0=x0\n", " X=x\n", " Pfinal=zeros(D,D)\n", " (gam, PP,xx,xx0, Nused, cond)=procedure_online_sampling(m,n,l,N,0,N1,alpha,tol,d,X0,X,C,0,0,zeros(n,n),Pfinal)\n", " PP1[Nphases]=PP\n", " xx0=xx0[:,1:Nused] \n", " xx=xx[:,1:Nused] \n", " L=cholesky(PP)\n", " L=L.U\n", " if d>1\n", " L=inversedlift(L,d,n)\n", " end\n", " Linv=inv(L)\n", " A1=[zeros(n,n) for i=1:m]\n", " for i=1:m\n", " A1[i]=A[i]\n", " end\n", " X0=x0\n", " X=x\n", " (u, uj, gam, xxx,xxx0,PP,lam,ubis)=proba_bounds(m,n,l,N1,beta,beta1,tol,d,X0,X,C)\n", " u1=zeros(1,phases)\n", " uj1=zeros(1,phases)\n", " gam1=zeros(1,phases)\n", " U1=zeros(1,phases)\n", " Nresult=zeros(1,phases+1)\n", " Nresult[1]=Nused\n", " condi1=zeros(2,1000,phases)\n", " while N1-Nused>convert(Int64, n*(n+1)/2)+N && Nphases1.05\n", " X=zeros(n,N1-Nused)\n", " X0=x0[:,1+Nused:N1]\n", " for j=1:N1-Nused\n", " b=rand((1:m),l)\n", " X[:,j]=Linv*X0[:,j] \n", " for k=1:l\n", " X[:,j]=A1[b[k]]*X[:,j]\n", " end\n", " X[:,j]=L*X[:,j] \n", " end\n", " for i=1:m\n", " A1[i]=L*A1[i]*Linv\n", " end \n", " i=Nphases\n", " ( gam1[i], PP1[i+1],xx,xx0,Nuse,cond1)=procedure_online_sampling(m,n,l,N,Nused,N1,alpha,tol,d, X0,X,C,xx0,xx,L,Pfinal)\n", " s=size(cond1)\n", " condi1[:,1:s[2],i]=cond1\n", " L=cholesky(PP1[i+1])\n", " Nused+=Nuse\n", " Nphases+=1\n", " Nresult[Nphases]=Nused\n", " L=L.U\n", " if d>1\n", " L=inversedlift(L,d,n)\n", " end\n", " Linv=inv(L)\n", " end\n", " if N1-Nused>0\n", " X=zeros(n,N1-Nused)\n", " X0=x0[:,1+Nused:N1]\n", " for j=1:N1-Nused\n", " b=rand((1:m),l)\n", " X[:,j]=Linv*X0[:,j] \n", " for k=1:l\n", " X[:,j]=A1[b[k]]*X[:,j]\n", " end\n", " X[:,j]=L*X[:,j] \n", " end\n", " for i=1:m\n", " A1[i]=L*A1[i]*Linv\n", " end \n", " i=Nphases\n", " (u1[i], gam1[i], PP1[i+1],U1[i],xx,xx0,ubis)=proba_bounds_offline_sampling(m,n,l,N1-Nused,beta,beta1,tol,d,X0,X,C,xx0,xx,Nused,L)\n", " Nphases+=1\n", " end\n", " return (u,u1[1:Nphases-1],U1[1:Nphases-1], sqrt(eigmax(PP)/eigmin(PP)),[sqrt(eigmax(PP1[i])/eigmin(PP1[i])) for i=1:Nphases],Nresult[1:Nphases],Nphases-1)\n", "end" ] }, { "cell_type": "code", "execution_count": 12, "id": "portuguese-exception", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "offline (generic function with 14 methods)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# function offline(Ntrials,N, N1, phases,A,m,n,l,d,beta, beta1, tol, 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 on the\n", "# JSR of a system A of m modes of dimension nxn, based on N randomly sampled trajectories of length l, \n", "# for the offline sampling with the ratio of samples used in phase 1, for several number of phases. \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 N1 contained\n", "# in the vector N1 and each value of phases, the average of the bounds over Ntrials runs is computed. The\n", "# upper bounds are based on [2, Thm. 5.1] and [2,Thm. 5.2] when d=1; on [2, Thm. 6.1] \n", "# otherwise.. \n", "# The true value of the JSR of the system is also computed.\n", "#\n", "# The default values of the arguments give rise to the first experiment of [2,Chapter 5.7.1]\n", "#\n", "# Inputs: -N: total number of samples\n", "# -N1: Array of positive integers containing all the ratios number of trajectories to be sampled in phase 1.\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", "# -phases: number of phases to use.\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(phases))x(size(N1)) float matrix. up[1,1] contains the average of the bound of [2,Thm. 3.1] over Ntrials runs, for N samples.\n", "# -up1: (size(phases))x(size(N1)) float matrix. up1[i,:] contains the average of the bound of [2,Thm. 5.1] over Ntrials runs, for each value of N1, for phases[i].\n", "# -UU: (size(phases))x(size(N1)) float matrix. U1[i,:] contains the average of the bound of [2,Thm. 5.2] over Ntrials runs, for each value of N1, for phases[i].\n", "# -P2: (size(phases))x(size(N1)) float matrix. P2[1,1] contains the average of the condition number of Psol based on N samples over Ntrials runs.\n", "# -P12: (size(phases))x(size(N1)) float matrix. P12[i,:] contains the average of the condition number of Pp over Ntrials runs, for each value of N1, for phases[i].\n", "\n", "\n", "function offline(Ntrials=10,N=2000, N1=[0.02,0.05,0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.85,0.9,0.95,0.98], phases=[1,2,3],A=generatesystem(2,8),m=2,n=8,l=1,d=1,beta=0.95, beta1=0.95, tol=0.001, C=1000)\n", "p=length(N1)\n", "D=length(phases)\n", "N1=N1*N\n", "N2=N*ones(p,1)-N1\n", "up1=zeros(D,p)\n", "up=zeros(D,p)\n", "P12=zeros(D,p)\n", "P2=zeros(D,p)\n", "gamma1=zeros(D+1,p)\n", "gamma=zeros(D,p)\n", "UU=zeros(D,p)\n", " uj1=0\n", "for i=1:Ntrials\n", " x0,x=generatetrajectories(m,n,l,A,N)\n", " for j=1:D\n", " for k=1:p\n", " X0=x0\n", " X=x\n", " (u, u1,U1,PP2,PP12)=algo_offline_sampling(m,n,l,convert(Int64,N1[k]),convert(Int64,floor(N2[k]/phases[j])),beta,beta1,tol,d,X0,X,A, C, phases[j])\n", " up1[j,k]+=u1[j]\n", " P12[j,k]+=PP12[j]\n", " UU[j,k]+=U1[j]\n", " if phases[j]==1\n", " X0=x0\n", " X=x\n", " (u,uj1)=algo_fast(m,n,l,convert(Int64,N1[k]),convert(Int64,N2[k]),beta,beta1,tol,d,X0,X, C) \n", " end\n", " end\n", " end\n", " X0=x0\n", " X=x\n", "(u, uj, gam,xxx,xxx0,PP,ubis)=proba_bounds(m,n,l,N,beta,beta1,tol,d,X0,X,C)\n", " up[1,1]+=u\n", " P2[1,1]+=sqrt(eigmax(PP)/eigmin(PP))\n", "end\n", "up1=up1/Ntrials\n", "up=up/Ntrials\n", "P12=P12/Ntrials\n", "P2=P2/Ntrials\n", "gamma1=gamma1/Ntrials\n", "gamma=gamma/Ntrials\n", "UU=UU/Ntrials\n", "jsr=trueJSR(A)\n", " N1=N1./N\n", " return (up,up1,P2,P12,UU,uj1)\n", "end" ] }, { "cell_type": "code", "execution_count": 13, "id": "herbal-tourist", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "online (generic function with 15 methods)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# function online(Ntrials,N1, N2, phases,A,alpha , m,n,l,d,beta, beta1, tol, 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 on the\n", "# JSR of a system A of m modes of dimension nxn, based on N2 randomly sampled trajectories of length l, \n", "# for the online sampling with N2. \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 N12 contained in the vector N2, the average of the bounds over Ntrials runs is computed. The\n", "# upper bounds are based on [2, Thm. 5.1] and [2,Thm. 5.2] when d=1; on [2, Thm. 6.1] \n", "# otherwise.. \n", "# The true value of the JSR of the system is also computed.\n", "#\n", "#\n", "# Inputs: -N1:increment of samples to use in each Procedure Online\n", "# -N2: Array of positive integers containing all the numbers 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", "# -phases: maximum number of phases to use.\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(N2))x(size(Ntrials)) float matrix. up[i,:] contains the bound of [2,Thm. 3.1] for N2[i] samples, at each trial.\n", "# -up1: (size(N2))x(size(Ntrials)) float matrix. up1[i,:] contains the bound of [2,Thm. 5.1] for N2[i] samples, at each trial.\n", "# -UU: (size(N2))x(size(Ntrials)) float matrix. U1[i,:] contains the bound of [2,Thm. 5.2] for N2[i] samples, at each trial.\n", "# -P2: (size(N2))x(size(Ntrials)) float matrix. P2[1,1] contain the condition number of Psol for N2[i] samples, at each trial.\n", "# -P12: (size(phases))x(size(N2))xsize(Ntrials) float matrix. P12[:,i,:] contains the condition number of Pp for N2[i] samples, at each trial, at each phase. \n", "# -Nres: (size(phases))x(size(N2))xsize(Ntrials) float matrix. Nres[:,i,:] contains the number of samples used in each phase for N2[i] samples, at each trial.\n", "# -Nphases: (size(N2))x(size(Ntrials)) float matrix. Number of phases used for N2[i] samples, at each trial..\n", "\n", "function online(Ntrials=10,N1=5, N2=100 .+300*[k for k=0:13],phases=10, A=generatesystem(2,8),alpha=0.1,m=2,n=8,l=1,d=1,beta=0.95, beta1=0.95, tol=0.001, C=1000)\n", "p=length(N2)\n", "up1=zeros(p,Ntrials)\n", "up=zeros(p,Ntrials)\n", "P12=zeros(phases+1,p,Ntrials)\n", "P2=zeros(p,Ntrials)\n", "res=zeros(phases,p,Ntrials)\n", "gamma1=zeros(p,Ntrials)\n", "gamma=zeros(p,Ntrials)\n", "Nphases=zeros(p,Ntrials)\n", "UU=zeros(phases,p,Ntrials)\n", "x0=zeros(n,N2[p])\n", "x=zeros(n,N2[p]) \n", "for i=1:Ntrials\n", " for k=1:p\n", " if k==1\n", " X0,X=generatetrajectories(m,n,l,A,N2[k])\n", " x0[:,1:N2[k]]=X0\n", " x[:,1:N2[k]]=X\n", " else\n", " X0,X=generatetrajectories(m,n,l,A,N2[k]-N2[k-1])\n", " x0[:,1+N2[k-1]:N2[k]]=X0\n", " x[:,1+N2[k-1]:N2[k]]=X\n", " end\n", " X0=x0\n", " X=x \n", " ( u, u1,U1,PP2,PP12,Nres,Npha)=algo_online_sampling(m,n,l,N1,N2[k],beta,beta1,tol,d,X0,X, C, A, phases,alpha)\n", " up1[k,i]=u1[end]\n", " up[k,i]=u\n", " P12[1:Npha+1,k,i]=PP12\n", " UU[1:Npha,k,i]=U1\n", " P2[k,i]=PP2\n", " Nphases[k,i]=Npha\n", " res[1:Npha,k,i]=Nres[1:Npha]\n", " end\n", "end\n", " return (up,up1,P2,P12,UU, res, Nphases)\n", "end\n" ] }, { "cell_type": "code", "execution_count": 14, "id": "behavioral-simulation", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "online_several_dimensions (generic function with 14 methods)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# function online_several_dimensions(Ntrials,N1, N2, phases,alpha,m,n,l,d,beta, beta1, tol, C)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function analyzes the probabilistic upper bounds on the\n", "# JSR of a system A of m modes of dimension nxn, based on N2 randomly sampled trajectories of length l, \n", "# for the online sampling, for several values of 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 in the vector n, the average of the bounds over Ntrials runs is computed. The\n", "# upper bounds are based on [2, Thm. 5.1] and [2,Thm. 5.2] when d=1; on [2, Thm. 6.1] \n", "# otherwise.\n", "# The true value of the JSR of the system is also computed.\n", "#\n", "#\n", "# Inputs: -N1:increment of samples to use in each Procedure Online\n", "# -N2: 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: Array of positive integers 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", "# -phases: maximum number of phases to use.\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(Ntrials)) float matrix. up[i,:] contains the bound of [2,Thm. 3.1] for N2[i] samples, at each trial.\n", "# -up1: (size(n))x(size(Ntrials)) float matrix. up1[i,:] contains the bound of [2,Thm. 5.1] for N2[i] samples, at each trial.\n", "# -UU: (size(n))x(size(Ntrials)) float matrix. U1[i,:] contains the bound of [2,Thm. 5.2] for N2[i] samples, at each trial.\n", "# -P2: (size(n))x(size(Ntrials)) float matrix. P2[1,1] contain the condition number of Psol for N2[i] samples, at each trial.\n", "# -P12: (size(phases))x(size(n))xsize(Ntrials) float matrix. P12[:,i,:] contains the condition number of Pp for N2[i] samples, at each trial, at each phase. \n", "# -Nres: (size(phases))x(size(n))xsize(Ntrials) float matrix. Nres[:,i,:] contains the number of samples used in each phase for N2[i] samples, at each trial.\n", "# -Nphases: (size(n))x(size(Ntrials)) float matrix. Number of phases used for N2[i] samples, at each trial..\n", "\n", "function online_several_dimensions(Ntrials=10,N1=5,N2=1000,phases=10,alpha=0.1,m=2,n=[3,4,5,6,7,8,10],l=1,d=1,beta=0.95, beta1=0.95, tol=0.001, C=1000)\n", "p=length(n)\n", "up1=zeros(phases,p,Ntrials)\n", "up=zeros(p,Ntrials)\n", "P12=zeros(phases+1,p,Ntrials)\n", "P2=zeros(p,Ntrials)\n", "pfinal=zeros(phases+1,p,Ntrials)\n", "gamma1=zeros(phases,p,Ntrials)\n", "gamma=zeros(p,Ntrials)\n", "UU=zeros(phases,p,Ntrials)\n", " Nphases=zeros(p,Ntrials)\n", " res=zeros(phases,p,Ntrials)\n", "for i=1:Ntrials\n", " for k=1:p\n", " D=binomial(d+n[k]-1,d)\n", " A=generatesystem(m,n[k])\n", " x0,x=generatetrajectories(m,n[k],l,A,N2)\n", " X0=x0\n", " X=x \n", " (u, u1,U1,PP2,PP12,Nres,Npha)=algo_online_sampling(m,n[k],l,N1,N2,beta,beta1,tol,d,X0,X, C, A, phases,alpha)\n", " up1[1:Npha,k,i]=u1\n", " up[k,i]=u\n", " P12[1:Npha+1,k,i]=PP12\n", " UU[1:Npha,k,i]=U1\n", " P2[k,i]=PP2\n", " Nphases[k,i]=Npha\n", " res[1:Npha,k,i]=Nres[1:Npha]\n", " end\n", "end\n", " return (up,up1,P2,P12,UU, res, Nphases)\n", "end\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 15, "id": "handled-israel", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "testheuristic (generic function with 1 method)" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ " #function heuristic(Ntrials,N,Nphainterpol,m,n,l,d,beta, beta1, tol, C)\n", "#\n", "# Author: Anne Rubbens (anne.rubbens@student.uclouvain.be)\n", "# Date: June 2021\n", "#\n", "# Description: This function analyzes the probabilistic upper 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", "# for the heuristic [2, Heuristic 1], for several values of 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 in the vector n, the average of the bounds over Ntrials runs is computed. The\n", "# upper bounds are based on [2, Thm. 5.1] and [2,Thm. 5.2] when d=1; on [2, Thm. 6.1] \n", "# otherwise.\n", "# The true value of the JSR of the system is also computed.\n", "#\n", "#\n", "# Inputs: -Nphainterpol:number of samples to be used in each phase i allowed.\n", "# -N: 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: Array of positive integers 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", "# -C: Positive integer describing the bound on the norm of P in Problem (1). \n", "#\n", "# Output: -up: (size(n))x(size(Ntrials)) float matrix. up[i,:] contains the bound of [2,Thm. 3.1] for N2[i] samples, at each trial.\n", "# -up1: (size(n))x(size(Ntrials)) float matrix. up1[i,:] contains the bound of [2,Thm. 5.1] for N2[i] samples, at each trial.\n", "# -UU: (size(n))x(size(Ntrials)) float matrix. U1[i,:] contains the bound of [2,Thm. 5.2] for N2[i] samples, at each trial.\n", "# -P2: (size(n))x(size(Ntrials)) float matrix. P2[1,1] contain the condition number of Psol for N2[i] samples, at each trial.\n", "# -P12: (size(phases))x(size(n))xsize(Ntrials) float matrix. P12[:,i,:] contains the condition number of Pp for N2[i] samples, at each trial, at each phase. \n", "function heuristic(Ntrials=10,N=2000, Nphainterpol=[45 -45\n", " 30 -150\n", " 20 -80\n", " 20 -100\n", " 10 -10\n", " 10 -60\n", " 20 -160],m=2,n=[3,4,5,6,7,8,10],l=1,d=1,beta=0.95, beta1=0.95, tol=0.001, C=1000)\n", "p=length(n)\n", "up1=zeros(10,p,Ntrials)\n", "up=zeros(p,Ntrials)\n", "upj1=zeros(10,p,Ntrials)\n", "upj=zeros(p,Ntrials)\n", "P1=zeros(10,p,Ntrials)\n", "P=zeros(p,Ntrials)\n", "P12=zeros(10,p,Ntrials)\n", "P2=zeros(p,Ntrials)\n", "pfinal=zeros(10,p,Ntrials)\n", "pfinalj=zeros(10,p,Ntrials)\n", "gamma1=zeros(10,p,Ntrials)\n", "gamma=zeros(p,Ntrials)\n", "UU=zeros(10,p,Ntrials)\n", " upbis1=zeros(10,p,Ntrials)\n", " upbis=zeros(p,Ntrials)\n", " jsr=zeros(p,1)\n", " Npha=zeros(100,1)\n", " phases=0\n", "\n", " for k=1:p\n", " D=binomial(d+n[k]-1,d)\n", " s=size(Nphainterpol)\n", " Npha[1:s[1]].=Nphainterpol[:,1]*D+Nphainterpol[:,2]\n", " if D <=5\n", " phases=5\n", " Npha[2:5]=D*(D+1)/2*ones(1,4)\n", " Np=N-sum([Npha[i] for i=1:4])\n", " p=5\n", " while Np0\n", " Npha[1]=N-2*(D*(D+1)/2+1)\n", " Npha[2]=D*(D+1)/2+1\n", " Npha[3]=D*(D+1)/2+1\n", " phases=3\n", " elseif N-(D*(D+1)/2+1)>0\n", " Npha[1]=N-(D*(D+1)/2+1)\n", " Npha[2]=D*(D+1)/2+1\n", " phases=2\n", " else\n", " Npha[1]=0.8*N\n", " Npha[2]=0.2*N\n", " phases=2\n", " end\n", " end\n", " elseif D<=10\n", " phases=6\n", " Np=N-sum([Npha[i] for i=1:5])\n", " p=6\n", " while Np0\n", " Npha[1]=N-2*(D*(D+1)/2+1)\n", " Npha[2]=D*(D+1)/2+1\n", " Npha[3]=D*(D+1)/2+1\n", " phases=3\n", " elseif N-(D*(D+1)/2+1)>0\n", " Npha[1]=N-(D*(D+1)/2+1)\n", " Npha[2]=D*(D+1)/2+1\n", " phases=2\n", " else\n", " Npha[1]=0.8*N\n", " Npha[2]=0.2*N\n", " phases=2\n", " end\n", " end\n", " else \n", " phases=10\n", " Npha[s[1]+1:phases]=ceil(D*(D+1)/4)*ones(1,phases-s[1])\n", " Np=N-sum([Npha[i] for i=1:9])\n", " p=10\n", " while Np0\n", " Npha[1]=N-(p-1)*(D*(D+1)/2+1)\n", " for i=2:p\n", " Npha[i]=D*(D+1)/2+1\n", " end\n", " phases=p\n", " else\n", " Npha[1]=0.8*N\n", " Npha[2]=0.2*N\n", " phases=2\n", " end\n", " end\n", " end\n", " A=generatesystem(m,n[k])\n", " jsr[k]=trueJSR(A)\n", " for i=1:Ntrials\n", " x0,x=generatetrajectories(m,n[k],l,A,N)\n", " X0=x0\n", " X=x \n", " (ubis,ubis1, u, u1,U1,PP2,PP12)=testheuristic(m,n[k],l,N,Npha,phases,beta,beta1,tol,d,x0,x, C, A)\n", " up1[1:phases-1,k,i]=u1\n", " upbis1[1:phases-1,k,i]=ubis1\n", " upbis[k,i]=ubis\n", " up[k,i]=u\n", " P12[1:phases,k,i]=PP12\n", " UU[1:phases-1,k,i]=U1\n", " P2[k,i]=PP2\n", " end\n", "end\n", " return (up,up1,P2,P12,UU,jsr,upbis,upbis1)\n", "end\n", "function testheuristic(m,n,l,N,Npha,phases,beta,beta1,tol,d,x0,x, C, A)\n", " D=binomial(d+n-1,d)\n", " PP1=[zeros(D,D) for i=1:phases]\n", " Nphases=1\n", " X0=x0\n", " X=x\n", " (u, uj, gam, xx,xx0,PP,U,ubis)=proba_bounds(m,n,l,convert(Int64,Npha[1]), beta,beta1,tol,d,X0,X,C)\n", " PP1[Nphases]=PP\n", " Nused=convert(Int64,Npha[1])\n", " xx0=xx0[:,1:Nused] \n", " xx=xx[:,1:Nused] \n", "L=cholesky(PP)\n", " L=L.U\n", " if d>1\n", " L=inversedlift(L,d,n)\n", " end\n", " Linv=inv(L)\n", " A1=[zeros(n,n) for i=1:m]\n", " for i=1:m\n", " A1[i]=A[i]\n", " end\n", " X0=x0\n", " X=x\n", " (u, uj, gam, xxx,xxx0,PP,Uinitulie, ubis)=proba_bounds(m,n,l,N,beta,beta1,tol,d,X0,X,C) \n", " u1=zeros(1,phases-1)\n", " uj1=zeros(1,phases-1)\n", " gam1=zeros(1,phases-1)\n", " U1=zeros(1,phases-1)\n", " upbis=zeros(1,phases-1)\n", " while Nphases1.05\n", " X=zeros(n,convert(Int64,Npha[Nphases+1]))\n", " X0=x0[:,1+Nused:convert(Int64,Npha[Nphases+1])+Nused]\n", " for j=1:convert(Int64,Npha[Nphases+1])\n", " b=rand((1:m),l)\n", " X[:,j]=Linv*X0[:,j] \n", " for k=1:l\n", " X[:,j]=A1[b[k]]*X[:,j]\n", " end\n", " X[:,j]=L*X[:,j] \n", " end\n", " for i=1:m\n", " A1[i]=L*A1[i]*Linv\n", " end \n", " i=Nphases\n", " (u1[i], gam1[i], PP1[i+1],U1[i],xx,xx0,upbis[i])=proba_bounds_offline_sampling(m,n,l,convert(Int64,Npha[Nphases+1]),beta,beta1,tol,d, X0,X,C,xx0,xx,Nused,L)\n", " L=cholesky(PP1[i+1])\n", " \n", " Nused+=convert(Int64,Npha[Nphases+1])\n", " Nphases+=1\n", " L=L.U\n", " if d>1\n", " L=inversedlift(L,d,n)\n", " end\n", " Linv=inv(L)\n", " end\n", " return (ubis,upbis,u,u1,U1, sqrt(eigmax(PP)/eigmin(PP)),[sqrt(eigmax(PP1[i])/eigmin(PP1[i])) for i=1:Nphases])\n", "end\n" ] } ], "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 }