From 1f12e47e6d1fe8b365806aa0b42de237970d53bf Mon Sep 17 00:00:00 2001 From: Nao Pross Date: Wed, 24 May 2023 16:50:06 +0200 Subject: Take deliverables for Robust MPC from yanzhen According to table 12 - generate_disturbances - simulate_uncertain - compute_tube_contorller - compute_minRPI - compute_tightening - MPC_TUBE - MPC_TUBE/eval (contained in MPC_TUBE.m) - MPC_TUBE_script (and its output MPC_TUBE_params.mat) --- templates/cprnd.m | 296 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 296 insertions(+) create mode 100755 templates/cprnd.m (limited to 'templates/cprnd.m') diff --git a/templates/cprnd.m b/templates/cprnd.m new file mode 100755 index 0000000..4183ac3 --- /dev/null +++ b/templates/cprnd.m @@ -0,0 +1,296 @@ +function [X S] = cprnd(N,A,b,options) +%CPRND Draw from the uniform distribution over a convex polytope. +% X = cprnd(N,A,b) Returns an N-by-P matrix of random vectors drawn +% from the uniform distribution over the interior of the polytope +% defined by Ax <= b. A is a M-by-P matrix of constraint equation +% coefficients. b is a M-by-1 vector of constraint equation +% constants. +% +% cprnd(N,A,b,options) allows various options to be specified. + +% 'method' Specifies the algorithm. One of the strings 'gibbs', +% 'hitandrun' (the default), and 'achr'. The default +% algorithm 'hitandrun' is a vanilla hit-and-run sampler +% [1]. 'gibbs' specifies a Gibbs sampler, 'achr' is the +% Adaptive Centering Hit-and-Run algorithm of [2]. +% +% 'x0' Vector x0 is a starting point which should be interior +% to the polytope. If x0 is not supplied CPRND uses the +% Chebyshev center as the initial point. +% +% 'isotropic' Perform an adaptive istropic transformation of the +% polytope. The values 0, 1 and 2 respectively turn +% off the transformation, construct it during a runup +% period only, or continuously update the +% tranformation throughout sample production. The +% transformation makes sense only for the Gibbs and +% Hit-and-run samplers (ACHR is invariant under +% linear transformations). +% +% 'discard' Number of initial samples (post-runup) to discard. +% +% 'runup' When the method is gibbs or hitandrun and the +% isotropic transformation is used, this is the +% number of initial iterations of the algorithm in +% the untransformed space. That is a sample of size +% runup is generated and its covariance used as the +% basis of a transformation. +% When the method is achr runup is the number of +% conventional hitandrun iterations. See [2]. +% +% 'ortho' Zero/one flag. If turned on direction vectors for +% the hitandrun algorithm are generated in random +% orthonormal blocks rather than one by +% one. Experimental and of dubious utility. +% +% 'quasi' Allows the user to specify a quasirandom number generator +% (such as 'halton' or 'sobol'). Experimental and of +% dubious utility. +% +% + +% By default CPRND employs the hit-and-run sampler which may +% exhibit marked serial correlation and very long convergence times. +% +% References +% [1] Kroese, D.P. and Taimre, T. and Botev, Z.I., "Handbook of Monte +% Carlo Methods" (2011), pp. 240-244. +% [2] Kaufman, David E. and Smith, Robert L., "Direction Choice for +% Accelerated Convergence in Hit-and-Run Sampling", Op. Res. 46, +% pp. 84-95. +% +% Copyright (c) 2011-2012 Tim J. Benham, School of Mathematics and Physics, +% University of Queensland. + + function y = stdize(z) + y = z/norm(z); + end + + p = size(A,2); % dimension + m = size(A,1); % num constraint ineqs + x0 = []; + runup = []; % runup necessary to method + discard = []; % num initial pts discarded + quasi = 0; + method = 'achr'; + orthogonal = 0; + isotropic = []; + + % gendir generates a random unit (direction) vector. + gendir = @() stdize(randn(p,1)); + + % Alternative function ogendir forces directions to come in + % orthogonal bunches. + Ucache = {}; + function u = ogendir() + if length(Ucache) == 0 + u = stdize(randn(p,1)); + m = null(u'); % orthonormal basis for nullspace + Ucache = mat2cell(m',ones(1,p-1)); + else + u = Ucache{end}'; + Ucache(end) = []; + end + end + + % Check input arguments + + if m < p+1 + % obv a prob here + error('cprnd:obvprob',['At least ',num2str(p+1),' inequalities ' ... + 'required']); + end + + if nargin == 4 + if isstruct(options) + + if isfield(options,'method') + method = lower(options.method); + switch method + case 'gibbs' + case 'hitandrun' + case 'achr' + otherwise + error('cprnd:badopt',... + ['The method option takes only the ' ... + 'values "gibbs", "hitandrun", and "ACHR".']); + end + end + + if isfield(options,'isotropic') + % Controls application of isotropic transformation, + % which seems to help a lot. + isotropic = options.isotropic; + end + + if isfield(options,'discard') + % num. of samples to discard + discard = options.discard; + end + + if isfield(options,'quasi') + % Use quasi random numbers, which doesn't seem to + % help much. + quasi = options.quasi; + if quasi && ~ischar(quasi), quasi='halton'; end + if ~strcmp(quasi,'none') + qstr = qrandstream(quasi,p,'Skip',1); + gendir = @() stdize(norminv(qrand(qstr,1),0,1)'); + end + end + + if isfield(options,'x0') + % initial interior point + x0 = options.x0; + end + + if isfield(options,'runup') + % number of points to generate before first output point + runup = options.runup; + end + + if isfield(options,'ortho') + % Generate direction vectors in orthogonal + % groups. Seems to help a little. + orthogonal = options.ortho; + end + + else + x0 = options; % legacy support + end + end + + % Default and apply options + + if isempty(isotropic) + if ~strcmp(method,'achr') + isotropic = 2; + else + isotropic = 0; + end + end + + if orthogonal + gendir = @() ogendir(); + end + + % Choose a starting point x0 if user did not provide one. + if isempty(x0) + x0 = chebycenter(A,b); % prob. if p==1? + end + + % Default the runup to something arbitrary. + if isempty(runup) + if strcmp(method,'achr') + runup = 10*(p+1); + elseif isotropic > 0 + runup = 10*p*(p+1); + else + runup = 0; + end + end + + % Default the discard to something arbitrary + if isempty(discard) + if strcmp(method,'achr') + discard = 25*(p+1); + else + discard = runup; + end + end + + X = zeros(N+runup+discard,p); + + n = 0; % num generated so far + x = x0; + + % Initialize variables for keeping track of sample mean, covariance + % and isotropic transform. + M = zeros(p,1); % Incremental mean. + S2 = zeros(p,p); % Incremental sum of + % outer prodcts. + S = eye(p); T = eye(p); W = A; + + while n < N+runup+discard + y = x; + + % compute approximate stochastic transformation + if isotropic>0 + if n == runup || (isotropic > 1 && n > runup) + T = chol(S,'lower'); + W = A*T; + end + y = T\y; + end + + switch method + + case 'gibbs' + % choose p new components + for i = 1:p + % Find points where the line with the (p-1) components x_i + % fixed intersects the bounding polytope. + e = circshift(eye(p,1),i-1); + z = W*e; + c = (b - W*y)./z; + tmin = max(c(z<0)); tmax = min(c(z>0)); + % choose a point on that line segment + y = y + (tmin+(tmax-tmin)*rand)*e; + end + + case 'hitandrun' + % choose a direction + u = gendir(); + % determine intersections of x + ut with the polytope + z = W*u; + c = (b - W*y)./z; + tmin = max(c(z<0)); tmax = min(c(z>0)); + % choose a point on that line segment + y = y + (tmin+(tmax-tmin)*rand)*u; + + case 'achr' + % test whether in runup or not + if n < runup + % same as hitandrun + u = gendir(); + else + % choose a previous point at random + v = X(randi(n),:)'; + % line sampling direction is from v to sample mean + u = (v-M)/norm(v-M); + end + % proceed as in hit and run + z = A*u; + c = (b - A*y)./z; + tmin = max(c(z<0)); tmax = min(c(z>0)); + % Choose a random point on that line segment + y = y + (tmin+(tmax-tmin)*rand)*u; + end + + if isotropic>0 + x = T*y; + else + x = y; + end + + X(n+1,:)=x'; + n = n + 1; + + % Incremental mean and covariance updates + delta0 = x - M; % delta new point wrt old mean + M = M + delta0/n; % sample mean + delta1 = x - M; % delta new point wrt new mean + if n > 1 + S2 = S2 + (n-1)/(n*n)*delta0*delta0'... + + delta1*delta1'; + S0 = S; + S = S2/(n-1); % sample covariance + else + S = eye(p); + end + + end + + X = X((discard+runup+1):(N+discard+runup),:); + +end -- cgit v1.2.1