Commit 8e1035e4 authored by Haoyue Yang's avatar Haoyue Yang
Browse files
parents 27c6e322 73524500
......@@ -32,7 +32,23 @@ classdef dd < export.Data
header = {};
plotfnc = @plot;
splitN = inf;
% If a function is not continuous, the jumping points of the
% function can be defined by this structure. Then, the export data
% will be manipulated so that, the jumps are plotted as jumps and
% not as ramps.
% The property must be a cell array with structures as entries.
% These structures must have the fields:
% *idx_t: The column index of the time vector in M
% *idx_j: The column index of the variable that has the jump in M
% *tj: the time where the jump occurrs.
% As an example: If the time vector has the index 1 and the
% function is saved in column 3 and the jumps occurrs at t=5, then
% the required structure is generated by:
% struct('idx_t', 1, 'idx_j', 3, 'tj', 5)
jumps = [];
addPoint = [];
end
methods
......@@ -54,6 +70,11 @@ classdef dd < export.Data
end
function o = out(obj)
if ~isempty(obj.addPoint)
addData = obj.M(obj.addPoint, :);
end
if obj.N > 0 && obj.N < size(obj.M, 1)
x = linspace(obj.M(1,1), obj.M(end,1), obj.N)';
o = interp1(obj.M(:,1), obj.M(:,2:end), x);
......@@ -80,6 +101,10 @@ classdef dd < export.Data
end
end
if ~isempty(obj.addPoint)
%% add the required points
o = sortrows([o; addData], 1);
end
end
function is = isempty(obj)
......
......@@ -46,8 +46,8 @@ classdef ddd < export.Data
end
function out = get.out(obj)
o.x = linspace(obj.xRange(1), obj.xRange(2), 60);
o.y = linspace(obj.yRange(1), obj.yRange(2), 60);
o.x = linspace(obj.xRange(1), obj.xRange(2), obj.N);
o.y = linspace(obj.yRange(1), obj.yRange(2), obj.N);
[o.x, o.y] = meshgrid(o.x, o.y);
[grid.x, grid.y] = meshgrid(obj.x, obj.y);
o.z = interp2(grid.x, grid.y, obj.z', o.x, o.y);
......
......@@ -24,12 +24,14 @@ function Psi = odeSolver(A, zdisc)
% - created on 03.08.2018 by Simon Kerschbaum
% preliminaries
RelTol = 1e-12;
AbsTol = 1e-12;
ndisc=length(zdisc);
if ~isa(A,'function_handle')
if size(A,1)~= ndisc
error('If passed as discretized function, then size(A,1) must equal length(zdisc)!')
end
opts = odeset('RelTol',1e-9,'AbsTol',1e-9,'vectorized','on');
opts = odeset('RelTol', RelTol,'AbsTol', AbsTol,'vectorized','on');
% if discretized version was passed, the solver can use the vectorized option!
if size(A,2)~= size(A,3)
error('Size(A,2) must equal size(A,3), such that A is quadratic!')
......@@ -50,12 +52,12 @@ if ~isa(A,'function_handle')
A_for_vec = @(z) squeeze(AInterp({z(:)}));
end
else
opts = odeset('RelTol',1e-9,'AbsTol',1e-9,'vectorized','off');
ADisc = misc.eval_pointwise(A,zdisc);
opts = odeset('RelTol', RelTol, 'AbsTol', AbsTol,'vectorized','off');
ADisc = misc.eval_pointwise(A,zdisc(1));
n = size(ADisc,2);
A_for_vec = @(z) misc.fDiag(A,n,z);
end
Psi_mu_vec(1:ndisc,1:ndisc,1:n^2) = 0; % preallocate
Psi_mu_vec(1:ndisc,1:ndisc,1:n^2) = nan; % preallocate
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
......
classdef Odes < handle & matlab.mixin.Copyable
%misc.ss.Odes (standing for ordinary differential equations) is class to
%implement multiple different state spaces (see doc/ss).
properties
ss (1, :) cell; % cell-array of state spaces
type (1, :) cell; % cell array of names of ss
C (1, :) cell; % cell of output gains / matrices
D (1, :) cell; % cell of feedthrough gains / matrices
end
properties (Dependent = true)
odes (1, 1) ss; % connection of all state spaces in the ss property.
% % The connection is performed by connecting InputName and
% % OutputName, hence specification of those properties of
% % all ss-cell-elements is essential.
numTypes (1, 1); % number of different types
OutputName (:, 1) cell; % OutputNames of ss - unique, sorted
InputName (:, 1) cell; % InputNames of ss - unique, sorted
n (1, :) double; % size(ss(:).A, 1)
end
methods
function obj = Odes(varargin)
%% misc.ss.Odes constructs an object containing ss
% Example:
% ss1 = ss(1, 1, 1, 1, 'InputName', 'u1', 'OutputName', 'y1');
% ss2 = ss(2, 2, 2, 2, 'InputName', 'u2', 'OutputName', 'y2');
% myOdes = misc.ss.Odes(ss1, 'ss1', ss2, 'ss2')
if nargin > 0
it = 1;
while it <= nargin
if isa(varargin{it}, 'ss') && ischar(varargin{it+1})
obj.type{end+1} = varargin{it+1};
[obj.ss{end+1}, obj.C{end+1}, obj.D{end+1}] ...
= obj.seperateOutputfromSs(varargin{it}, obj.type{end});
it = it + 2;
else
error('input must be name value pairs of (ss, char-array)');
end
end
obj.verifyProperties();
end % if nargin > 0
end % misc.ss.Odes() Constructor
function verifyProperties(obj)
% verification of properties of misc.ss.Odes. Checks length of ss and
% type and that type is unique.
assert(misc.isunique(obj.type), 'type of misc.Odes must be unique');
assert(isequal(size(obj.ss), size(obj.type)), ...
'number of ss and types are not equal');
assert(isequal(size(obj.ss), size(obj.C)), ...
'number of ss and C are not equal');
assert(isequal(size(obj.ss), size(obj.D)), ...
'number of ss and D are not equal');
end % verifyProperties()
function obj = add(obj, varargin)
% with add() multiple misc.ss.Odes are combined
% Example:
% ss1 = ss(1, 1, 1, 1, 'InputName', 'u1', 'OutputName', 'y1');
% ss2 = ss(2, 2, 2, 2, 'InputName', 'u2', 'OutputName', 'y2');
% Ode1 = misc.ss.Odes(ss1, 'ss1');
% Ode2 = misc.ss.Odes(ss2, 'ss2');
% Odes12 = Ode1.add(Ode2);
for myOde = varargin{:}
assert(isa(myOde, 'misc.ss.Odes'), 'input must be an misc.ss.Odes-object');
for it = 1 : myOde.numTypes
assert(~any(contains(obj.type, myOde.type{it})), ...
['The ss-type ', myOde.type{it}, ' is already defined']);
obj.ss{end+1} = myOde.ss{it};
obj.C{end+1} = myOde.C{it};
obj.D{end+1} = myOde.D{it};
obj.type{end+1} = myOde.type{it};
end
end
obj.verifyProperties();
end % add()
function obj = exchange(obj, myType, newSs)
% Exchange one ss with another one
% Example:
% ss1 = ss(1, 1, 1, 1, 'InputName', 'u1', 'OutputName', 'y1');
% ss2 = ss(2, 2, 2, 2, 'InputName', 'u2', 'OutputName', 'y2');
% Ode = misc.ss.Odes(ss1, 'ss1');
% Ode.exchange('ss1', ss2)
[obj.ss{obj.index(myType)}, obj.C{obj.index(myType)}, obj.D{obj.index(myType)}] ...
= obj.seperateOutputfromSs(newSs, myType);
obj.verifyProperties();
end % exchange()
function c = plus(a, b)
% Example:
% ss1 = ss(1, 1, 1, 1, 'InputName', 'u1', 'OutputName', 'y1');
% ss2 = ss(2, 2, 2, 2, 'InputName', 'u2', 'OutputName', 'y2');
% Ode1 = misc.ss.Odes(ss1, 'ss1');
% Ode2 = misc.ss.Odes(ss2, 'ss2');
% Odes12 = Ode1 + Ode2
if isempty(a)
c = b;
elseif isempty(b)
c = a;
elseif isa(a, 'misc.ss.Odes') && isa(b, 'misc.ss.Odes')
c = copy(a);
c = c.add(b);
else
error('input data type not recognized');
end
a.verifyProperties();
end %plus()
function obj = remove(obj, myType)
% remove specified ss and type
% Example:
% ss1 = ss(1, 1, 1, 1, 'InputName', 'u1', 'OutputName', 'y1');
% ss2 = ss(2, 2, 2, 2, 'InputName', 'u2', 'OutputName', 'y2');
% myOdes = misc.ss.Odes(ss1, 'ss1', ss2, 'ss2');
% myOdes.remove('ss1')
if ~iscell(myType)
myType = {myType};
end
for it = 1 : numel(myType)
selectRemaining = (1:1:obj.numTypes) ~= obj.index(myType{it});
obj.ss = obj.ss(selectRemaining);
obj.C = obj.C(selectRemaining);
obj.D = obj.D(selectRemaining);
obj.type = obj.type(selectRemaining);
end
obj.verifyProperties();
end % remove()
function idx = index(obj, myType)
% find ss by naming its type and return index-value
% Example:
% ss1 = ss(1, 1, 1, 1, 'InputName', 'u1', 'OutputName', 'y1');
% ss2 = ss(2, 2, 2, 2, 'InputName', 'u2', 'OutputName', 'y2');
% myOdes = misc.ss.Odes(ss1, 'ss1', ss2, 'ss2');
% myOdes.index('ss2')
idx = [];
for it = 1 : obj.numTypes
if isequal(obj.type{it}, myType)
idx = it;
end
end
end % index()
function thisSs = getSs(obj, myType)
% get ss specified by type
% Example:
% ss1 = ss(1, 1, 1, 1, 'InputName', 'u1', 'OutputName', 'y1');
% ss2 = ss(2, 2, 2, 2, 'InputName', 'u2', 'OutputName', 'y2');
% myOdes = misc.ss.Odes(ss1, 'ss1', ss2, 'ss2');
% myOdes.getSs('ss2');
thisSs = obj.ss{obj.index(myType)};
end % getSs()
function thisC = getC(obj, myType)
% get ss specified by type
% Example:
% ss1 = ss(1, 1, 1, 3, 'InputName', 'u1', 'OutputName', 'y1');
% ss2 = ss(2, 2, 2, 4, 'InputName', 'u2', 'OutputName', 'y2');
% myOdes = misc.ss.Odes(ss1, 'ss1', ss2, 'ss2');
% myOdes.getC('ss2');
thisC = obj.C{obj.index(myType)}.D; % yes, .D and not .C is correct
end % getC()
function thisN = getN(obj, myType)
% get ss specified by type
% Example:
% ss1 = ss(1, 1, 1, 1, 'InputName', 'u1', 'OutputName', 'y1');
% ss2 = ss(2, 2, 2, 2, 'InputName', 'u2', 'OutputName', 'y2');
% myOdes = misc.ss.Odes(ss1, 'ss1', ss2, 'ss2');
% myOdes.getN('ss2');
thisN = size(obj.ss{obj.index(myType)}.A, 1);
end % getN()
function thisD = getD(obj, myType)
% get ss specified by type
% Example:
% ss1 = ss(1, 1, 1, 1, 'InputName', 'u1', 'OutputName', 'y1');
% ss2 = ss(2, 2, 2, 2, 'InputName', 'u2', 'OutputName', 'y2');
% myOdes = misc.ss.Odes(ss1, 'ss1', ss2, 'ss2');
% myOdes.getD('ss2');
thisD = obj.D{obj.index(myType)}.D;
end % getD()
function result = isempty(obj)
% A misc.ss.Odes is assumed to be empty if numTypes is zero or empty
% Example:
% isempty(misc.ss.Odes())
result = isempty(obj.numTypes) || (obj.numTypes == 0);
end %isempty()
function result = isequal(A, B, varargin)
% isequal compares all parameters of A and B and varargin and returns
% true if they are equal, and false if not. Two misc.ss.Odes are assumed
% to be equal, if all ss and type are equal.
result = (A.numTypes == B.numTypes) && isequal(A.type, B.type);
if result
for it = 1 : A.numTypes
if result
result = isequal(A.ss(it), B.ss(it));
end
end
end
if result && (nargin > 2)
result = result && B.isequal(varargin{:});
end
end % isequal()
function [ic, stateName] = getInitialCondition(obj, varargin)
% getInitialCondition reads varargin if there are initial conditions
% defined as name-value-pair with the name according to obj.stateName and
% returns 2 cell arrays:
% 1. ic is a cell array of the values of the states initial condition
% 2. stateName is a cell array of the names of those states.
icStruct = misc.struct(varargin{:});
% Create initial condition for all states specified by obj.type
% If a value is specified in icStruct, than this value is used.
% Otherwise, a default value is considered.
icTemp = cell(obj.numTypes, 1);
for it = 1 : obj.numTypes
if misc.isfield(icStruct, obj.type{it})
icTemp{it} = misc.getfield(icStruct, obj.type{it});
else
icTemp{it} = zeros(size(obj.ss{it}.A, 1), 1);
end
end
ic = [icTemp(:)];
stateName = [obj.type(:)];
end % getInitialCondition()
%% get methods for Dependent properties
function myOutputName = get.OutputName(obj)
myData = [];
for it = 1 : obj.numTypes
myData = vertcat(myData, obj.ss{it}.OutputName, ...
obj.C{it}.OutputName, obj.D{it}.OutputName);
end
myOutputName = unique(myData);
end % get.OutputName();
function myInputName = get.InputName(obj)
myData = [];
for it = 1 : obj.numTypes
myData = vertcat(myData, obj.ss{it}.InputName, obj.D{it}.InputName);
end
myInputName = unique(myData);
end % get.InputName();
function myOdes = get.odes(obj)
if isempty(obj)
myOdes = [];
else
myOdes = misc.ss.connect(obj.InputName, obj.OutputName, ...
obj.ss{:}, misc.ss.parallel(obj.C{:}, obj.D{:}));
end % if
end % get.odes();
function numTypes = get.numTypes(obj)
numTypes = numel(obj.ss);
end %get.numTypes()
function n = get.n(obj)
n = zeros(size(obj.type));
for it = 1 : numel(n)
n(it) = size(obj.ss{1}.A, 1);
end % for it = 1 : numel(n)
end %get.n()
end % methods
methods (Static = true)
function [mySs, C, D] = seperateOutputfromSs(mySs, type)
% separateCfromSs() removes the output matrix of a state space model mySs
% and replaces it with the identity matrix, in order to obtain the states
% of the ode as output variables. The old output matrix C results in a
% separate state space
if all(strcmp(misc.ss.removeEnumeration(mySs.OutputName), type))
assert(misc.iseye(mySs.C) && isequal(size(mySs.C), size(mySs.A)), ...
['output of states space model does not include all states']);
C = ss();
D = ss();
else
% get C as state space
C = ss([], [], [], mySs.C, 'OutputName', mySs.OutputName);
C = misc.ss.setSignalName(C, 'input', {type}, size(C, 2));
% get D as state space
if ~isempty(mySs.D) && any(mySs.D(:) ~= 0)
D = ss([], [], [], mySs.D, ...
'OutputName', mySs.OutputName, 'InputName', mySs.InputName);
else
D = ss();
end
% change C and D in mySs
mySs = set(mySs, 'C', eye(size(mySs.A)), 'D', zeros(size(mySs.A, 1), size(mySs.B, 2)));
mySs = misc.ss.setSignalName(mySs, 'output', {type}, size(mySs.A, 1));
end % if
end % function [ss, C] = seperateCfromSs(ss)
end % methods (Static = true)
methods (Access = protected)
% Override copyElement method:
function cpObj = copyElement(obj)
% Make a shallow copy of all properties
cpObj = copyElement@matlab.mixin.Copyable(obj);
end
end % methods (Access = protected)
end
......@@ -20,7 +20,7 @@ end
% clean state-spaces
for it = 1 : numel(inputSs)
% remove feedthrouhgs
% remove feedthroughs
inputSs{it} = misc.ss.removeInputOfOutput(ss(inputSs{it}));
inputSs{it}.OutputName = misc.ss.removeSingularEnumeration(inputSs{it}.OutputName);
inputSs{it}.InputName = misc.ss.removeSingularEnumeration(inputSs{it}.InputName);
......
function W = gramian(sys, t, t0)
% GRAMIAN controllability matrix.
% W = gramian(sys, t, t0) computes the gramian controllability matrix of
% the state space system sys. For each timestep t. If t0 is not set, the
% default value is 0. For the computation the result from
%
% C. Van Loan, "Computing integrals involving the matrix exponential," in IEEE Transactions on Automatic Control, vol. 23, no. 3, pp. 395-404, June 1978.
% doi: 10.1109/TAC.1978.1101743
% URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1101743&isnumber=241661
%
% are used. In particular, the formula below (1.4) is applied. The gramian
% is given as:
%
% W(t0, t1) = int_t0^t1 expm(A*(tau-t0)) * b * b^T * expm(A^T(tau-t0)) dtau
% = int_t0^t1 expm(A*(t1 - tau) * b * b^T * expm(A^T(t1 - tau)) dtau
%
if nargin == 2
t0 = 0;
end
n = size(sys.A, 1);
Q = [-sys.A sys.B*sys.B'; ...
zeros(n) sys.A'] * quantity.Discrete( t - t0, 'grid', t, 'gridName', 't');
temp = expm(Q);
F3 = temp(n+1:end,n+1:end);
G2 = temp(1:n,n+1:end);
W = F3' * G2;
end
......@@ -31,7 +31,7 @@ s12 = parallel(s1, s2, inIndex{1}, inIndex{2}, outIndex{1}, outIndex{2});
% recursive call for parallel connection of further state space models in varargin
if ~isempty(varargin)
s12 = parallel(s12, varargin{:});
s12 = misc.ss.parallel(s12, varargin{:});
end
end
\ No newline at end of file
function [u, y, x] = planTrajectory(sys, t, varargin)
% PLANTRAJECTORY Computes a trajectory for a lumped-parameter systems
% [u, y, x] = planTrajectory(sys, t, varargin) computes a transition
% u(t) : x(t0) = x0 --> x(t1) = x1
% for a state-space system sys of the form
% d/dt x = A x + bu.
% Initial and end value x0 and x1 can be defined by name value pairs.
% Default value is zero. The computation of the transition is based on
% the Gramian controllability matrix
% W(t0, t1) = int_t0^t1 Phi(t0, tau) b b^T Phi^T(t0, tau) d tau.
% The used formula is
% u(t) = -b^T Phi^T(t0, t) W^(-1)(t0, t1) ( x0 - Phi(t0, t1) x1 )
% x(t) = Phi(t0, t) * W(t0, t) * W0^(-1)(t0, t1) * ( Phi(t0, t1) * x1 - x0)
% y(t) = C x(t)
%
p = misc.Parser();
p.addParameter('x0', zeros(size(sys.A, 1), 1));
p.addParameter('x1', zeros(size(sys.A, 1), 1));
p.parse(varargin{:});
x0 = p.Results.x0;
x1 = p.Results.x1;
% prepare time vectors
t = t(:);
t0 = t(1);
t1 = t(end);
%% Compute the state transition matrix: Phi(t,tau) = expm(A*(t-tau) )
% Att0 = sys.A * quantity.Symbolic( sym('t') - sym('tau'), 'grid', {t, t}, 'gridName', {'t', 'tau'} );
Phi_t0 = expm(sys.A * quantity.Discrete( t - t0, 'grid', {t}, 'gridName', {'t'} ));
invPhi_t1 = expm(sys.A * quantity.Discrete( t1 - t, 'grid', {t}, 'gridName', {'t'} ));
%% compute the gramian controllability matrix
W1 = misc.ss.gramian(sys, t, t0);
W1_t1 = W1.at(t1);
% Formel aus dem Bernstein:
u = sys.b.' * invPhi_t1.' / W1_t1 * (x1 - Phi_t0.at(t1) * x0);
% Berechnung der
x = Phi_t0 * x0 + W1 * invPhi_t1' / W1_t1 * (x1 - Phi_t0.at(t1) * x0);
%
y = sys.c * x;
%% Alternative Lösung nach Chen:
% W0 = int( invPhi_t0 * sys.B * sys.B' * invPhi_t0', 't', t0, 't');
% W0 = invPhi_t0 * W1 * invPhi_t0';
% u1 = -sys.B' * invPhi_t0' / W0.at(t1) * (x0 - Phi_t1.at(t0) * x1);
% xW0 = Phi_t0 * ( x0 - W0 / W0.at(t1) * (x0 - Phi_t1.at(t0) * x1));
%
% % simulation results have shown, does the average of both soultions leads to a better than each.
% u = (u0 + u1) / 2;
......@@ -75,6 +75,10 @@ classdef Gain < handle & matlab.mixin.Copyable
end
end % plus()
function myGains = add(a, b)
myGains = a + b;
end % add()
function obj = exchange(obj, newGain)
if all(strcmp(obj.inputType, newGain.inputType))
for it = 1 : numel(newGain.outputType)
......@@ -150,10 +154,16 @@ classdef Gain < handle & matlab.mixin.Copyable
if ~iscell(outputNames)
outputNames = {outputNames};
end
outputNames = outputNames(...
cellfun(@(v) any(strcmp(v, obj.outputType.')), outputNames));
if isempty(outputNames)
newObj = misc.Gain();
else
newObj = misc.Gain(obj.inputType, ...
obj.valueOfOutput(outputNames), ...
'outputType', outputNames, ...
'lengthOutput', obj.lengthOfOutput(outputNames));
end
end % gainOfOutput()
function thisValue = valueOfOutput(obj, outputNames)
......@@ -172,6 +182,28 @@ classdef Gain < handle & matlab.mixin.Copyable
end
end % valueOfOutput(obj, outputNames)
function obj = removeOutput(obj, outputType)
myOutputRows = getOutputRows(obj, outputType);
if ~isempty(myOutputRows)
myOutputRowsSelector = true(sum(obj.lengthOutput), 1);
myOutputRowsSelector(myOutputRows) = false;
newLengthOutput = obj.lengthOutput(~strcmp(obj.outputType, outputType));
newOutputType = obj.outputType(~strcmp(obj.outputType, outputType));
newValue = obj.value(myOutputRowsSelector, :);
if isempty(newOutputType)
obj = misc.Gain();
else
obj = misc.Gain(obj.inputType, newValue, ...
'outputType', newOutputType, ...
'lengthOutput', newLengthOutput);
end
% obj.value = obj.value(myOutputRowsSelector, :);
% obj.outputType = obj.outputType(~strcmp(obj.outputType, outputType));
% obj.lengthOutput = lengthOutputBackup(~strcmp(obj.outputType, outputType));
obj.verifySizes();
end
end % removeOutput()
function thisLength = lengthOfOutput(obj, varargin)
if nargin > 1
thisOutputType = varargin{1};
......@@ -193,6 +225,22 @@ classdef Gain < handle & matlab.mixin.Copyable
assert(isequal([sum(obj.lengthOutput), obj.lengthInput], size(obj.value)));
end % verifySizes(obj)
function obj = strrepOutputType(obj, oldText, newText)
% replace strings in type
assert(ischar(oldText) || isstring(oldText), 'oldText must be a char-array or string');
assert(ischar(newText) || isstring(newText), 'newText must be a char-array or string');
for it = 1 : numel(obj.outputType)
obj.outputType{it} = strrep(obj.outputType{it}, oldText, newText);
end
end % strrepOutputType()
function obj = strrepInputType(obj, oldText, newText)
% replace strings in type
assert(ischar(oldText) || isstring(oldText), 'oldText must be a char-array or string');
assert(ischar(newText) || isstring(newText), 'newText must be a char-array or string');
obj.inputType = strrep(obj.inputType, oldText, newText);
end % strrepInputType()
function OutputName = get.OutputName(obj)
OutputName = cell(sum(obj.lengthOutput), 1);
myCounter = 1;
......@@ -225,8 +273,12 @@ classdef Gain < handle & matlab.mixin.Copyable
function result = isequal(A, B, varargin)
% isequal compares all parameters of A and B and varargin and returns
% true if they are equal, and false if not.
if (numel(A) == 0) || (numel(B) == 0)
result = (numel(A)==numel(B));
else
result = isequal(A.value, B.value) & strcmp(A.inputType, B.inputType) ...
& all(strcmp(A.outputType, B.outputType));
end
if result && (nargin > 2)
result = isa(B, 'misc.Gain') & B.isequal(varargin{:});
end
......