search

BobAubouin / Robust-and-Resilient-Control

This repository host the code developed during my master thesis: "Toward Robust and Resilient Cyber-Physical Systems: State awareness and Control under Physical faults and Cyber-attacks" at the IMS laboratory
MIT License
11 stars 1 forks source link

Hello Bob, #1

Closed klix0 closed 2 years ago

klix0 commented 2 years ago

my name is Klinto and I am trying to implement a new way to calculate mRPI sets in Matlab. I have seen that you have worked with this topic and I would like to ask you if you can help me with your knowledge on Matlab. Please replay me if you have time.

Have a nice day!

BobAubouin commented 2 years ago

Hello Klinto,

Yes I have used MPT and CORA toolboxe to compute zonotopic mRPI sets with Rakovic method. What's your issue ?

klix0 commented 2 years ago

Hello Bob,

thank you for your response! I am trying to implement the calculation of mRPI according to Schulze using MPT3 in Matlab too. Can you give me a personal email address or another way to talk, so I can explain you me work in more details. Of course if you have time.

Best Klinto

BobAubouin commented 2 years ago

Hello,

You can contact me on my academic mail : bob.aubouin-pairault@gipsa-lab.fr . I don't know if I will have the time to solve your issue but I can have a look !

Best, Bob

klix0 commented 2 years ago

Hello,

Thank you very much for your email!

Of course if you have time, just take a look. I have implemented the calculation of the Rakovic method in Matlab using the MPT3 toolbox and it works perfectly. Now I am trying to implement another method from Moritz Schulze, which is more effizient, but I have difficulties writing it in Matlab. I have written it in Julia and now I have to convert, but I have had some errors. I will share with you here the paper and the code in Julia.

If you are able to fix the problem, can you do it for me please? If it is necessary, I can pay for your work and for the time you will spend.

Best Klinto

using LazySets using LinearAlgebra using ControlSystems using Plots using InvariantSets using MathematicalSystems using Polyhedra

add https://github.com/ueliwechsler/InvariantSets.jl

NOTE: is defined in the package InvariantSets.ρ_matrix see (

https://github.com/ueliwechsler/InvariantSets.jl/blob/680cb87bd61cb26f3c1cfe97eaf27009d38ac016/src/helper/lazy_sets.jl#L65 )

otherwise the code for the logic should be self contained.

Functions needed

function explicit_maximal_RPI_set(Hd, bd, 𝒲, S, k) Hr0 = Hd br0 = bd Hrk = Array{Array}(undef, k+1) brk = Array{Vector}(undef, k+1) supk = Array{Vector}(undef, k+1) supk[1] = zeros(length(br0)) Hrk[1] = Hr0 brk[1] = br0 for i=2:k+1 supk[i] = supk[i-1] + InvariantSets.ρ_matrix(HdS^(i-2), 𝒲) Hrk[i] = [Hrk[i-1]; HdS^k] brk[i] = [brk[i-1]; bd - supk[i]] end return Hrk[end], brk[end] end """ 𝒟 = (𝕏,𝕌) = {x∈ℝⁿ | Hdx≤bd}""" function s_min_sufficient( S, 𝒲, Hd, bd; ε=0.1, verbose=false, k_max=1000 ) Hw, bw = tosimplehrep(𝒲) k = 0 while true k += 1 cond11 = (1 + ε)InvariantSets.ρ_matrix(HwS^k, 𝒲) .<= εbw bc = (1 + ε)sum(InvariantSets.ρ_matrix(HdS^j, 𝒲) for j in 0:k-1) cond12 = bc .<= bd (verbose) && @show k, cond11, cond12 if all(cond11) && all(cond12) return k, bc end k == k_max && error("Convergence not achieved for s = $k_max\n" "Check if the origin ∈ interior(𝒲)") end end function minimal_RPI_set_schulze(S, 𝒲, 𝕏::LazySet, 𝕌::LazySet, K; ε=0.1, k_max=1000, remove_redundancy=true, avoid_equality_check=false, verbose=true ) Hx, bx = tosimplehrep(𝕏) Hu, bu = tosimplehrep(𝕌) Hd = [Hx; Hu*K] bd = [bx; bu] @time s, bc = s_min_sufficient(S, 𝒲, Hd, bd; ε=ε, k_max=k_max) @info "Darup-Schulze: maximal s=$s" Hc = Hd 𝒞 = HPolytope(Hc, bc)

maximal_RPI_set does converge in at most s-1 steps

@time Hrk, brk = explicit_maximal_RPI_set(Hc, bc, 𝒲, S, s) 𝒫∞ = HPolytope(Hrk, brk) @show s return 𝒫∞, 𝒞, s end

Example Ueli

A = [1.0 1.0; 0.0 1.0] B = [0.5; 1.0] n = size(A,1) m = size(B,2) Q = I R = I K = dlqr(A,B,Q,R) S = A - B*K Hx = [Matrix{Float64}(I,n,n); -Matrix{Float64}(I,n,n)] bx = [4.0; 2.0; 8.0; 4.0] 𝒳 = HPolytope(Hx, bx) Hu = [Matrix{Float64}(I,m,m); -Matrix{Float64}(I,m,m)] bu = [1.0; 1.0] 𝒰 = HPolytope(Hu, bu) #or Hyperrectangle? 𝒲 = BallInf(zeros(n), 0.1) ε1 = 0.001 ε_mrpi1 = 0.1 Pinfty1, 𝒞1, k1 = minimal_RPI_set_schulze(S, 𝒲, 𝒳, 𝒰, K; ε=ε1); plot(𝒳) plot!(𝒞1, color=:gray,alpha=1) plot!(Pinfty1, color=:yellow,alpha=1)

=========================================================

Example from Anjeda Gaga:

=========================================================

note that comments use hashtags not %

make your own discrete linear system (plant)

A = [1 1; 0 1] # matrices in julia do not have a comma for row separation B = hcat([1; 1]) sys = @.** x⁺ = Ax + B*u

constraints on state Xc and input Uc

Xc_vertex = [2.0 -2.0; 2.0 2.0; -10.0 2.0; -10.0 -2.0]

Uc_vertex = hcat([1.0; -1.0]) # hcat is needed to create a matrix, since

VPolytope needs a matrix as input

Xc = VPolytope(Xc_vertex')

Xc = BallInf(zeros(2), 3.0) Uc = BallInf(zeros(1), 3.0)

Uc = VPolytope(Uc_vertex)

construct a convex set of system noise (2dim here)

W_vertex = [0.25 0.25; 0.25 -0.25; -0.25 -0.25; -0.25 0.25]

W = VPolytope(W_vertex)

set boundary of system noise which corresponds to W.

w_max = 0.25 W = BallInf(zeros(2), 0.25) plot(W) Q = diagm([1.0, 1.0]) R = 0.1 K = -ControlSystems.dlqr(A,B,Q,R)

closed loop system

Ak = A + B*K ε1 = 0.001 Pinfty1, 𝒞1, k1 = minimal_RPI_set_schulze(Ak, W, Xc, Uc, -K; ε=ε1) # note the minus for -K plot(Xc) plot!(𝒞1, color=:gray,alpha=1) plot!(Pinfty1, color=:yellow,alpha=1) [image: image.png]

On Thu, Jun 23, 2022 at 9:01 AM Bob Aubouin--Pairault < @.***> wrote:

Hello,

You can contact me on my academic mail : @.*** . I don't know if I will have the time to solve your issue but I can have a look !

Best, Bob

— Reply to this email directly, view it on GitHub https://github.com/Bobepsa/Robust-and-Resilient-Control/issues/1#issuecomment-1164029556, or unsubscribe https://github.com/notifications/unsubscribe-auth/AYPBACS5JST6SKK4BXGJT5TVQQDUZANCNFSM5ZNCKTAA . You are receiving this because you authored the thread.Message ID: @.***>

clear close all clc

% % Example 3 in Table 1 Schulze-Darup, Schaich & Cannon:

A = [.5,0; .5,1]; B = [1;0]; K = dlqr(A,B,eye(2),1e-2); A = A-B*K; D = -eye(2);

Lambda = [-K; K; eye(2); -eye(2)]; lambda = [.2;.3;.7;.5;.3;.5];

G = [speye(2);-speye(2)]/.05;

P = A'A; V = [1,-1,1,-1; 1,1,-1,-1]/20; cand = zeros(1,4); for i = 1:4 cand(i) = V(:,i)'P*V(:,i); end R = sqrt(max(cand));

decayRate = sqrt(1-max(eig(A))); lambdaMin = min(svd(A)); Fact = R/(1-decayRate)sqrt(1/lambdaMin); cand = zeros(1,length(lambda)); for i = 1:length(lambda); cand(i) = Factsqrt(Lambda(i,:)PLambda(i,:)'); end hTilde = max(cand);

LambdaLifted = [Lambda,zeros(size(Lambda,1),1); zeros(2),[-1;1/1.5]]; lambdaLifted = [lambda;0;1];

[TermA,TermB,~] = MRPIFullDimDist(LambdaLifted,lambdaLifted,G,A,D,200,'MPC book Example');

Num = 10; InvSetA = cell(1,Num); InvSetb = cell(1,Num); iters = zeros(1,Num);

alpha = linspace(0,1.5,Num); for i = 1:Num [InvSetA{i},InvSetb{i},iters(i)] = MRPIfixedAlpha(Lambda,lambda,A,D,G,200,alpha(i)); end

figure(1) hold on for i = 1:Num plot(Polyhedron(InvSetA{i},InvSetb{i}),'alpha',.2) end

x0 = [.25;-.5]; x = [x0,zeros(2,50)]; for i = 1:50 x(:,i+1) = A*x(:,i); end

plot(x(1,:),x(2,:),'bx') hold off

BobAubouin commented 2 years ago

Hey, Can you send me the code and the paper by mail please ? I'm not used to Julia, I can only have look on matlab but please send me all the necessary files (I believe for instance that MRPIFullDimDist function is not define in the code you send me).

Best, Bob