%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2023, Amon Lahr, Simon Muntwiler, Antoine Leeman & Fabian Flürenbrock Institute for Dynamic Systems and Control, ETH Zurich.
%
% All rights reserved.
%
% Please see the LICENSE file that has been included as part of this package.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [H_tube,h_tube,n_iter] = compute_minRPI(K_tube,params)
    % YOUR CODE HERE
    A = params.model.A+params.model.B*K_tube;
    % A = params.model.A;
%     display(A)
    nx = params.model.nx;
    omega = Polyhedron('A',[eye(nx);-eye(nx)],'b',zeros(1,2*nx)');
    % display(omega.A)
    % display(omega.b)
    n_iter = 0;
    Hw = params.constraints.DisturbanceMatrix;
    hw = params.constraints.DisturbanceRHS;
    W = Polyhedron('A',Hw,'b',hw); 
    while true 
        omega_last = copy(omega);
        % display(n_iter)
        % A_last = omega.A;
        % b_last = omega.b;
        display(n_iter);
        display(omega.b)
        omega = omega.plus(A^(n_iter)*W);
        
        if eq(omega.minHRep(),omega_last.minHRep())
        % if eq(omega,omega_last)
            break
        end
        n_iter = n_iter+1;
    end
    H_tube = omega.A;
    h_tube = omega.b;
end