Commit 964b1f59 authored by Patrick Göttsch's avatar Patrick Göttsch
Browse files

Gitlab Runner Update from lnsi-auxfiles/master:Update .gitlab-ci.yml

parent eb447dfb
% Exercise 7.5 :: Least Squares Identification
%
% TUHH :: Institut for Control Systems :: Control Systems Theory and Design
% Last update: 23.01.2009
clear all
clc
load cs7_LSsysdat
Ts = 1; % sampling time
n = 3; % model order
M = mkm(y3,u3,n);
yn = y3(n+1:length(y3)); % FIXME
p = pinv(M'*M)*M'*yn; % get parameters % FIXME
den = [1 p(1:n)']; % FIXME
num = p(n+1:2*n)'; % FIXME
G = tf(num,den,Ts);
[p z] =pzmap(G);
%% Plot simulation results of this model with the other signals
Tfin = length(u3)-1;
T = (0:Ts:Tfin)';
figure(1) % step
[y,T,x]=lsim(G,u1,T);
subplot(211), plot(T,[y y1]), title('step input : true plant and model'), legend('model','plant')
subplot(212), plot(T,y-y1), title('step input : error between true plant and Model')
figure(2) % sin
[y,T,x]=lsim(G,u2,T);
subplot(211), plot(T,[y y2]), title('sinusoidal input : true plant and model'), legend('model','plant')
subplot(212), plot(T,y-y2), title('sinusoidal input : error between true plant and Model ')
figure(3) % white noise
[y,T,x]=lsim(G,u3,T);
subplot(211), plot(T,[y y3]), title('white noise input : true plant und Model'), legend('model','plant')
subplot(212), plot(T,y-y3), title('white noise input : error between true plant and model')
% Exercise 7.5 :: Least Squares Identification - inverses
%
% TUHH :: Institut for Control Systems :: Control Systems Theory and Design
% Last update: 23.01.2009
clear all
clc
format short e
load cs7_LSsysdat
for i=1:10
M = mkm(y2,u2,i);
svd(M'*M)
end
for i=1:20
M = mkm(y3,u3,i);
svd(M'*M)
end
function M = mkm(y,u,n)
% M = mkm(y,u,n)
% Form the measurement matrix M
% y - measured outputs
% u - measured inputs
% n - order of the system
ny = length(y);
nu = length(u);
ndata = min(nu,ny)-1;
for i=n:ndata
M(i-n+1,1:n) = -y(i:-1:i-n+1)';
M(i-n+1,n+1:2*n) = u(i:-1:i-n+1)';
end
In Problem 3.1 there is a mismatch between the weights that were calculated in subtask a) by hand and thoose computed by Matlab in subtask b).
Like some of you already have observed correctly this is because of a preprocessing step that is performed by Matlab automatically.
Both datasets (input data and target data) are normalized to the interval [-1,1] in order to improve the quality of the trained network.
Attached is a Matlab script that compares the original output data and that obained by a NN with the weights obtained in b).
The preprocessing which was performed by Matlab is not considered. One can observe that output is some sort of a scaled version of the original data.
To receive exatly the same weights that were computed in a) with Matlab, one has to export the Network to the workspace e.g. as 'net'. After applying the lines
>> net.inputs{1}.processFcns{1} = 'removeconstantrows';
>> net.outputs{2}.processFcns{1} = 'removeconstantrows';
it can be reimported and trained. Please compare the Matlab Documentation.
https://de.mathworks.com/help/deeplearning/ug/choose-neural-network-input-output-processing-functions.html
\ No newline at end of file
% LNSI Exercise 3_1: Manually hand tuning of the weights
% Exercise class at :http://webcast.tuhh.de/Mediasite/Play/e6e3e15dee944bf79e183bdf47c175931d
clear
clc
%% Load data
load('approx.mat')
%% Plot original curve
figure()
plot(phi,g_phi,'DisplayName','Original Data','LineWidth',2)
grid on
legend
hold on
%% Grid input-axis
x=-10:0.1:10;
%% Weight choice 1
% CONDITIONS ON WEIGHTS derived in class
%w_1_10=7*w_1_11;
%w_1_11<0
%w_1_20=-5*w_1_21;
w_2_10=-1;
w_2_11=2;
w_2_12=8;
w_1_11=-3; % governs slope at the first inflection point.
w_1_21=0.5;% governs slope at the second inflection point
w_1_10=7*w_1_11;
w_1_20=-5*w_1_21;
% compute output and plot
y1=w_2_11./(1+ exp(-1*(w_1_10+w_1_11*x))) + w_2_12./(1+ exp(-(w_1_20+w_1_21*x)))+w_2_10;
figure()
plot(phi,g_phi,'DisplayName','Original Data','LineWidth',2)
grid on
legend
hold on
plot(x,y1,'--','DisplayName','Weight Choice 1','LineWidth',1)
%% Weight choice 2
% CONDITIONS ON WEIGHTS derived in class
%w_1_10=7*w_1_11;
%w_1_11<0
%w_1_20=-5*w_1_21;
w_2_10=1;
w_2_11=-2;
w_2_12=8;
w_1_11=3; % governs slope at the first inflection point.
w_1_21=0.5; % governs slope at the second inflection point.
w_1_10=7*w_1_11;
w_1_20=-5*w_1_21;
% compute output and plot
y2=w_2_11./(1+ exp(-1*(w_1_10+w_1_11*x))) + w_2_12./(1+ exp(-(w_1_20+w_1_21*x)))+w_2_10;
figure()
plot(phi,g_phi,'DisplayName','Original Data','LineWidth',2)
grid on
legend
hold on
plot(x,y2,'-.','DisplayName','Weight choice 2','LineWidth',1.5)
%% Trained network Weights taken from nntoolbox for visualization
% Input to Layer1
w_1_11=0.5;
w_1_21=-3;
% Bias to Layer1
w_1_10=-2.5;
w_1_20=-21;
% Hidden Layer to Output layer
w_2_11=2.1811;
w_2_12= 0.54528;
% Bias to Output layer
w_2_10=-1.0157;
% Compute and plot
y3=w_2_11./(1+ exp(-1*(w_1_10+w_1_11*x))) + w_2_12./(1+ exp(-(w_1_20+w_1_21*x)))+w_2_10;
figure()
plot(phi,g_phi,'DisplayName','Original Data','LineWidth',2)
grid on
legend
hold on
plot(x,y3,'DisplayName','Scaled by nn tool','LineWidth',1)
% ph_init
if exist('ysim') clear ysim;
end
if exist('usim') clear usim;
end
stepint = 7; % initial setpoint
TankModel_Initialisation; % Initialize pH plant
Tsim = Tosam;
%% Simulate Open Loop Response to Excitation Signal
sim('TankModel_openloop',Tsim);
% Take half of the samples for training
N = floor(length(ysim.time)/2);
t_train = ysim.time(1:N);
u_train = usim.signals.values(1:N);
y_train = ysim.signals.values(1:N,:);
% Take the other half of the samples for validation
t_validate = ysim.time(N:2*N);
u_validate = usim.signals.values(N:2*N);
y_validate = ysim.signals.values(N:2*N,:);
% Plot data
figure;
subplot(211)
plot(t_validate-t_validate(1),y_validate)
ylabel('y [pH]')
xlabel('t [s]')
grid
subplot(212)
plot(t_validate-t_validate(1),u_validate)
ylabel('u [m^3/s]')
xlabel('t [s]')
grid
figure;
hist(y_train)
%% Investigate Time Constants
% Generate several open loop simulations to find out about the range of
% time constants in the nonlinear model. The definition of the sampling
% time is based on this
Tosam = 100; % number of sample
tsam = 0.20; % sampling time
stepint = 7; % initial setpoint
dsteptime = 500; % time of step disturbance
TankModel_Initialisation; % Initialize pH plant
%% Simulate Open Loop Response
figure;
% umin = 0.0;
% umax = 1.0;
% umin_mdl = 0.0;
% umax_mdl = 1.0;
% scale_range = 1.0;
Tsim = 50;
t1 = 0:tsam:Tosam-tsam;
uMag1 = 0;
uMag2 = 1;
uMag = [uMag1, uMag2];
%g = uMagnitude*ones(1,Tosam/tsam); % Excitation signal
% Setup figure
figure(1);
ylabel('y [pH]')
xlabel('t [s]')
grid
hold on
simmodel = 'TankModel_openloop';
% sim('TankModel_openloop_train',Tsim);
%% Plot step responses for initial values
g1 = uMag1*ones(1,Tosam/tsam); % Excitation signal
simin = [t1', g1'];
sim(simmodel,Tsim);
% Take simulation data and plot it
N = length(ysim.time);
t = ysim.time(1:N);
u = usim.signals.values(1:N);
y = ysim.signals.values(1:N,:);
yend = y(end);
% Plot data
figure(1);
plot(t-t(1),y,'b'); hold on
g2 = uMag2*ones(1,Tosam/tsam); % Excitation signal
simin = [t1', g2'];
sim(simmodel,Tsim);
% Take simulation data and plot it
N = length(ysim.time);
t = ysim.time(1:N);
u = usim.signals.values(1:N);
y = ysim.signals.values(1:N,:);
yend = [yend, y(end)];
% Plot data
figure(1);
plot(t-t(1),y,'b'); hold on
%% Bisection
ii = 1;
while(abs(yend(ii,1) - yend(ii,2)) > 0.1)
uMagMid = (uMag1 + uMag2)/2;
gk = uMagMid*ones(1,Tosam/tsam); % Excitation signal
simin = [t1', gk'];
sim(simmodel,Tsim);
% Take simulation data and plot it
N = length(ysim.time);
t = ysim.time(1:N);
u = usim.signals.values(1:N);
y = ysim.signals.values(1:N,:);
% Plot data
figure(1);
plot(t-t(1),y,'b'); hold on
if(y(end) > 7)
uMag1 = uMagMid;
%uMag2 = uMag2;
yend(ii+1,:) = [y(end), yend(ii,2)];
else
%uMag1 = uMag1;
uMag2 = uMagMid;
yend(ii+1,:) = [yend(ii,1), y(end)];
end
uMag = [uMag; [uMag1, uMag2]];
ii = ii+1;
end
%% Investigate Time Constants
% Generate several open loop simulations to find out about the range of
% time constants in the nonlinear model. The definition of the sampling
% time is based on this
Tosam = 100; % number of sample
tsam = 0.20; % sampling time
stepint = 7; % initial setpoint
dsteptime = 500; % time of step disturbance
TankModel_Initialisation; % Initialize pH plant
%% Simulate Open Loop Response
figure;
% umin = 0.0;
% umax = 1.0;
% umin_mdl = 0.0;
% umax_mdl = 1.0;
% scale_range = 1.0;
Tsim = 50;
t1 = 0:tsam:Tosam-tsam;
uMagnitude = cell(3,1);
% uMagnitude{1,1} = [-0.1:0.01:0.07]'; % Random magnitude of excitation signal
% uMagnitude{2,1} = [0.071:0.001:0.079]';
% uMagnitude{3,1} = [0.07999:1e-6:0.08001]';
% uMagnitude{4,1} = [0.081:0.001:0.089]';
% uMagnitude{5,1} = [0.09:0.010:1.000]';
% colors = {'b','g','r','g','b'};
% uMagnitude{1,1} = [0:0.05:0.45]'; % Random magnitude of excitation signal
% uMagnitude{2,1} = [0.480:0.002:0.499]';
% uMagnitude{3,1} = [0.5015:0.0001:0.5039]';
% uMagnitude{4,1} = [0.504:0.002:0.518]';
% uMagnitude{5,1} = [0.52:0.050:1.000]';
% colors = {'b','g','r','g','b'};
uMagnitude{1,1} = [0:0.01:0.46]'; % Random magnitude of excitation signal
uMagnitude{2,1} = [0.460:0.0025:0.530]';
uMagnitude{3,1} = [0.530:0.010:1.000]';
colors = {'b','r','b'};
%g = uMagnitude*ones(1,Tosam/tsam); % Excitation signal
% Setup figure
figure(1);
ylabel('y [pH]')
xlabel('t [s]')
grid
hold on
uydata = [];
for ii = 1:size(uMagnitude,1)
for jj = 1:size(uMagnitude{ii,1})
gk = uMagnitude{ii,1}(jj)*ones(1,Tosam/tsam); % Excitation signal
simin = [t1', gk'];
%sim('TankModel_openloop_train',Tsim);
sim('TankModel_openloop',Tsim);
% Take simulation data and plot it
N = length(ysim.time);
t = ysim.time(1:N);
u = usim.signals.values(1:N);
y = ysim.signals.values(1:N,:);
uydata = [uydata; [u(end), y(end)]];
% Plot data
figure(1);
plot(t-t(1),y,colors{ii}); hold on
end
end
figure(1);
xlabel('Time');
ylabel('pH Value');
figure(2);
plot(uydata(:,1),uydata(:,2));
xlabel('Reagent flow u');
ylabel('pH Value');
%% Multi-Sine Excitation Signal
% In the 'SINE' case, the sinusoids are chosen from the frequency grid
% freq = 2*pi*[1:Grid_Skip:fix(P/2)]/P intersected with pi*[BAND(1) BAND(2)].
% (for Grid_Skip see below.) For multi-input signals, the different inputs
% use different frequencies from this grid. An integer number of
% full periods is always delivered. The selected frequencies are obtained
% as [U,FREQS] = IDINPUT(....), where row ku of FREQS contains the
% the frequencies of input number ku. The resulting signal is affected by a
% 5th input argument SINEDATA:
% U = IDINPUT(N,TYPE,BAND,LEVELS,SINEDATA)
% where
% SINEDATA= [No_of_Sinusoids, No_of_Trials, Grid_Skip],
% meaning that No_of_Sinusoids are equally spread over the indicated
% BAND, trying No_of_Trials different, random, relative phases,
% until the lowest amplitude signal is found.
% Default SINEDATA = [10,10,1];
Tosam = 10000; % number of sample
tsam = 0.25; % sampling time
stepint = 7; % initial setpoint
TankModel_Initialisation; % Initialize pH plant
% Generate Multi-Sine Signal
gamma0 = 0.1;
gamma1 = 0.2;
gamma = 0.2;
% First Signal Portion - High Amplitude, Low Frequency
bandwidth = [1/(timcon)*0; 1/(timcon)*0.1]; % 0 to 0.3 times the system's open-loop time constant
N = 15; % Number of sinusoids
randomPhases = 10; % Number of random phases
levels = ([-1, 1])*gamma0 + [0.5, 0.5];
sinedata = [N, randomPhases, 1];
g0 = idinput(Tosam/tsam, 'Sine', bandwidth, levels, sinedata)';
% First Signal Portion - High Amplitude, Low Frequency
bandwidth = [1/(timcon)*0.1; 1/(timcon)*0.3]; % 0 to 0.3 times the system's open-loop time constant
N = 15; % Number of sinusoids
randomPhases = 10; % Number of random phases
levels = ([-1, 1])*gamma1;
sinedata = [N, randomPhases, 1];
g1 = idinput(Tosam/tsam, 'Sine', bandwidth, levels, sinedata)';
% Second Signal Portion - Low Amplitude, High Frequency
bandwidth = [1/(timcon)*0.3; 1/(2*tsam)]; % 0.5 to 3 times the system's open-loop time constant
N = 10; % Number of sinusoids
levels = [-1, +1]*gamma;
sinedata = [N, randomPhases, 1];
g2 = gamma*idinput(Tosam/tsam, 'Sine', bandwidth, levels, sinedata)';
gk = g0+g1+g2;
gk(gk < 0) = 0;
gk(gk > 1) = 1;
t1 = 0:tsam:Tosam-tsam;
simin = [t1', gk'];
plot(gk); % show input data for training
\ No newline at end of file
function x = MultilevelPRBS(simt,samt,timcon,alpha,levec)
% generate N-samples-constant signal
% simt : Time to simulate
% samt : Sampling time
% alpha : probability of u(t-1)
nu = round(simt/samt);
xinput = randn(1,nu);
% sim('randomsim',simt);
% xinput = rx'+meax;
nr = 0;
ticon = timcon/samt;
d = 1;
xoutput = [];
r = length(levec);
d = floor(r*rand(1)+1);
xoutput(1) = levec(d);
for i = 2:nu,
k = rand(1);
if (k > alpha)&(nr > ticon),
d = floor(r*rand(1)+1);
xoutput(i) = levec(d);
nr = 0;
else
xoutput(i) = xoutput(i-1);
nr = nr+1;
end
end
x = xoutput;
% -------------------------------> APCINIT.M <------------------------------
% ---------- Switches -----------
regty ='apc'; % Controller type (apc, pid, none)
refty ='myref'; % Reference signal (siggener/<var. name>)
simul ='nnet'; % Control object spec. (simulink/matlab/nnet)
if exist('simulink')~=5,
simul ='matlab'; % Simulink not present
end
% ---------- Initializations -----------
Ts = 0.2; % Sampling period (in seconds)
samples = 300 ; % Number of samples in simulation
u_0 = 0; % Initial control input
y_0 = 0; % Initial output
ulim_min = -Inf; % Minimum control input
ulim_max = Inf; % Maximum control input
% -- System to be Controlled (SIMULINK) --
integrator= 'ode45'; % Name of dif. eq. solver (f. ex. ode45 or ode15s)
sim_model = 'spm1'; % Name of SIMULINK model
% --- System to be Controlled (MATLAB) --
mat_model = 'springm'; % Name of MATLAB model
model_out = 'smout'; % Output equation (function of the states)
x0 = [0;0]; % Initial states
% ----- Neural Network Specification ------
% "nnfile" must contain the following variables which together define
% a NNARX model:
% NN, NetDef, W1, W2
% (i.e. regressor structure, architecture definition, and weight matrices)
nnfile = 'forward2'; % Name of file containing neural network model
% ------------ Reference filter ---------------
Am = [1]; % Denominator of filter
Bm = [1]; % Numerator of filter
% ---------- APC initializations -----------
N1 = 1; % Minimum prediction horizon (typically=nk)
N2 = 7; % Maximum prediction horizon (>= nb)
Nu = 2; % Control horizon
rho = 0.03; % Weight factor on differenced control signal
% ----------- Reference signal -------------
dc = 0; % DC-level
sq_amp = 3; % Amplitude of square signals (row vector)
sq_freq = 0.05; % Frequencies of square signals (column vector)
sin_amp = [0]; % Amplitude of sine signals (row vector)
sin_freq= [0]'; % Frequencies of sine signals (column vector)
Nvar = 0'; % Variance of white noise signal
myref =[0.3*ones(50,1);-0.3*ones(50,1);];
myref =[myref;1*ones(50,1);-1*ones(50,1)];
myref=[myref;2*ones(50,1);-2*ones(50,1)];
% -- Constant Gain Controller Parameters --
K=0.8; % PID parameters
Td=0.8; % PID
alf=0.1; % PID
Wi=0.2; % PID (1/Ti)