jens/HendiControlFirmware
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

- updated

Browse Source 
git-svn-id: http://moon:8086/svn/projects/HendiControl@5 fda53097-d464-4ada-af97-ba876c37ca34
This commit is contained in:
jens
2018-09-10 16:51:37 +00:00
parent c4e48007f1
commit 8fca4cd8c2
1 changed files with 134 additions and 50 deletions
Download Patch File Download Diff File Expand all files Collapse all files
matlab/brew_eval.m
+134 -50
View File
@@ -18,14 +18,13 @@
function brew_eval(varargin)
%
% brew_eval('M', 20.0, 'cont_ctrl', 0, 'kp', 0.1, 'ki', 0.1, 'kd', 100, 'T', 1800, 'mass_kleak', 0.2, 'P_q', 100, 'mass_delay', 100, [50 60 70 80]);
% brew_eval('M', 20.0, 'cont_ctrl', 1, 'kp', 2.0, 'ki', 0.02, 'kd', 400, 'T', 1800, 'mass_kleak', 0.5, 'P_q', 100, 'mass_delay', 60, 'dt', 1.00, [30 50 70]);
WITH_LOWPASS = 0;
WITH_RUNNING_MEAN = 0;
WITH_RATE_CONTROLLER = 1;
WITH_HYST = 0.0;
params = struct ('cont_ctrl',0, 'C',4.19e3, 'M',10, 'mass_kleak',0.5, 'mass_delay',50, 'P_min',0, 'P_max',3000, 'P_q',100, 'T',100, 'dt',1.0, 'Td',10.0, 'theta_lock',0.5, 'theta_amb',20.0, 'kp',10, 'ki',0.04, 'kd',500, 'pid_kleak',0.1, 'k_noise', 0.0, 'k_noise2', 0.0, 'f_noise2', 0.01);
units = struct ('cont_ctrl', '', 'C', 'J/(kg*K)', 'M', 'kg', 'mass_kleak', '', 'mass_delay', 's', 'P_min', 'W', 'P_max', 'W', 'P_q','W', 'T', 's', 'dt', 's', 'Td', 's', 'theta_lock', '°C', 'theta_amb', '°C', 'kp', '', 'ki', '', 'kd', '', 'pid_kleak', '', 'k_noise', '', 'k_noise2', '', 'f_noise2', 'Hz');
params = struct ('cont_ctrl',0, 'C',4.19e3, 'M',10, 'mass_kleak',0.5, 'mass_delay',50, 'P_min', 0, 'P_max',3000, 'P_q',100, 'T',100, 'dt',1.0, 'Td',10.0, 'theta_lock', 0.5, 'theta_amb',20.0, 'kp',10, 'ki',0.04, 'kd',500, 'pid_kleak',0.001, 'k_noise', 0.0, 'heatrate', 1.0);
units = struct ('cont_ctrl', '', 'C', 'J/(kg*K)', 'M', 'kg', 'mass_kleak', '', 'mass_delay', 's', 'P_min', 'W', 'P_max', 'W', 'P_q','W', 'T', 's', 'dt', 's', 'Td', 's', 'theta_lock', '°C', 'theta_amb', '°C', 'kp', '', 'ki', '', 'kd', '', 'pid_kleak', '', 'k_noise', '', 'heatrate', '°C/min');
% Parse parameters
names = fieldnames(params);
@@ -48,45 +47,121 @@ theta_last = params.theta_amb;
thetas = varargin{nargin}
heatRate = 0;
heatRate_f = 0;
heatRateMeasureInterval = 10.0; % seconds
heatRateMeasureInterval = 1.0; % seconds
heatRateMeasureCount = 0; % seconds
n = 0;
nSteps = params.P_max / params.P_q;
N = 2*fix(length(thetas)*params.T/params.dt);
maz = 17;
lp_N = 17;
[lp_b, lp_a] = butter (lp_N-1, 0.125);
maz = 7;
lp_N = 5;
[lp_b, lp_a] = butter (lp_N-1, 0.25);
lp_z = zeros(1, length(lp_b)-1);
lp_rate_N = 17;
[lp_rate_b, lp_rate_a] = butter (lp_rate_N-1, 0.125);
lp_rate_z = zeros(1, length(lp_rate_b)-1);
P_cont = 0;
y_rate = 0;
y = 0;
stbl_rate_z = [];
stbl_hold_z = [];
ovenState = 'HOLD';
controllerState = 'TRK';
curr = struct('kp', 0, 'ki', 0, 'kd', 0);
curr_rate = struct('kp', 0, 'ki', 0, 'kd', 0);
err0_hold = 0;
err0_rate = 0;
a_err = 0.5;
P_max = 0;
heatRate_soll = 0;
for k=1:length(thetas)
theta_soll = thetas(k);
hold_theta = 0;
temp_ok = 0;
P_max = 0;
printf('%d: Target temperatur %d °C\n', k, theta_soll);
while n < N,
err = theta_ist - theta_soll;
[y, yc, s_pid] = jpid(s_pid, params.kp, params.ki*params.dt, params.kd*params.dt, 1-params.pid_kleak, -err);
s_pid.y_max = 50;
P_rate = params.P_max;
% Rate controller
err_rate = heatRate - 1.0;
[y_rate, yc_rate, s_pid_rate] = jpid(s_pid_rate, 0.4*params.dt, 0.001*params.dt, 8.0*params.dt, 0.999, -err_rate);
s_pid_rate.y_max = 4000;
if WITH_RATE_CONTROLLER
if (hold_theta == 0)
P_rate = min(params.P_max, params.P_max*y_rate);
end
err1 = theta_ist - theta_soll;
err1_rate = heatRate - heatRate_soll;
err0_hold = (1-a_err)*err0_hold + a_err*err1;
err0_rate = (1-a_err)*err0_rate + a_err*err1_rate;
[err_rate_stbl, stbl_rate_z] = isStable(stbl_rate_z, err0_rate, 0.05);
[err_hold_stbl, stbl_hold_z] = isStable(stbl_hold_z, err0_hold, 0.05);
ovenState_n = ovenState;
controllerState_n = controllerState;
switch ovenState
case {'HEAT'}
heatRate_soll = params.heatrate;
curr_rate.kp = 4.0;
curr_rate.ki = 0.02;
curr_rate.kd = 4.0;
P_max = 0.5*P_max + 0.5*P_cont;
switch controllerState
case {'ACQ'}
if err_rate_stbl && (abs(err0_rate) < 0.1)
controllerState_n = 'TRK';
end
case {'TRK'}
curr_rate.kp = 8.0;
curr_rate.ki = 0.03;
curr_rate.kd = 8.0;
if abs(err0_rate) > 0.5
controllerState_n = 'ACQ';
end
endswitch
if abs(err0_hold) < 1.0
ovenState_n = 'HOLD';
controllerState_n = 'ACQ';
end
case {'HOLD'}
heatRate_soll = 0;
curr.kp = params.kp;
curr.ki = params.ki;
curr.kd = params.kd;
switch controllerState
case {'ACQ'}
if err_hold_stbl && (abs(err0_hold) < 0.1)
controllerState_n = 'TRK';
end
case {'TRK'}
curr.kp = params.kp;
curr.ki = params.ki;
curr.kd = params.kd;
if abs(err0_hold) > 0.2
controllerState_n = 'ACQ';
end
endswitch
if err0_hold < -1.0
if WITH_RATE_CONTROLLER
ovenState_n = 'HEAT';
end
controllerState_n = 'ACQ';
end
endswitch
if strcmp(controllerState, controllerState_n) == 0 || strcmp(ovenState, ovenState_n) == 0
printf('%d: State %s::%s -> %s::%s\n', n, ovenState, controllerState, ovenState_n, controllerState_n);
end
% Rate controller
ovenState = ovenState_n;
controllerState = controllerState_n;
P_cont = max(params.P_min, min(P_rate, P_rate*y));
P_cont_q = P_cont; %fix(P_norm * nSteps)*params.P_q;
% Hold controller
if strcmp(ovenState, 'HEAT')
[y, yc, s_pid] = jpid(s_pid, curr_rate.kp, curr_rate.ki*params.dt, curr_rate.kd/params.dt, 1-params.pid_kleak, -err0_rate);
P_cont = max(params.P_min, min(params.P_max, params.P_max*y));
s_pid.y_max = 80;
else
[y, yc, s_pid] = jpid(s_pid, curr.kp, curr.ki*params.dt, curr.kd/params.dt, 1-params.pid_kleak, -err0_hold);
P_cont = max(params.P_min, min(P_max, params.P_max*y));
s_pid.y_max = 80;
end
P_cont_q = fix(P_cont/params.P_q + 0.5)*params.P_q;
if params.cont_ctrl
P = P_cont_q;
else
@@ -95,43 +170,34 @@ for k=1:length(thetas)
P = params.P_max * heat;
end
[theta, s_mass] = mass(s_mass, params.dt, params.C, params.M, params.mass_kleak, P, params.mass_delay, theta_ist-params.theta_amb);
theta_ist = params.theta_amb + theta + params.k_noise*randn()/sqrt(12) + params.k_noise2*cos(params.f_noise2*2*pi*n*params.dt);
theta_ist = params.theta_amb + theta + params.k_noise*randn()/sqrt(12);
if WITH_RUNNING_MEAN
[theta_ist, maz] = movingavg(theta_ist, maz);
end
if WITH_LOWPASS
[theta_ist, lp_z] = filter(lp_b, lp_a, theta_ist, lp_z);
end
if heatRateMeasureCount <= 0
heatRate = 60*(theta_ist-theta_last)/heatRateMeasureInterval;
heatRateMeasureCount = heatRateMeasureInterval;
theta_last = theta_ist ;
else
heatRateMeasureCount = heatRateMeasureCount - params.dt;
theta_last = theta_ist;
end
heatRateMeasureCount = heatRateMeasureCount - params.dt;
[heatRate_f, lp_rate_z] = filter(lp_rate_b, lp_rate_a, heatRate, lp_rate_z);
n = n + 1;
theta_ist = theta_ist;
heatRate_(n) = heatRate_f;
theta_(n) = theta_ist;
err_(n) = err;
err_(n) = err0_hold;
p_(n) = P;
pc_(n) = P_cont;
pcq_(n) = P_cont_q;
y_(n) = y;
P_rate_(n) = P_rate;
err_rate_(n) = err_rate;
y_rate_(n) = y_rate;
err_rate_(n) = err0_rate;
if hold_theta > 0
hold_theta = hold_theta - 1;
if hold_theta == 0
disp('Temperature next')
break;
end
elseif abs(err) < params.theta_lock
elseif abs(err0_hold) < params.theta_lock
hold_theta = params.T/params.dt;
temp_ok = 1;
disp('Temperature is OK!')
end
end
end
@@ -141,17 +207,17 @@ t = (0:n-1)*params.dt;
subplot(3,1,1)
plot(t, theta_); grid; xlabel('t/s'); ylabel('Theta/°'); title('Temperature');
subplot(3,1,2)
plot(t, p_, 'b-', t, y_, 'g-', pc_, 'r-'); grid; xlabel('t/s'); ylabel('P/W'); title('Power Control'); legend('P','y','P_{cont}');
plot(t, p_, 'b-', t, pcq_, 'g-', t, pc_, 'r-'); grid; xlabel('t/s'); ylabel('P/W'); title('Power Control'); legend('P','y','P_{cont}');
subplot(3,1,3)
plot(t, err_); grid; xlabel('t/s'); ylabel('err/K'); title('Error'); legend('Error');
plot(t, err_, t, err_rate_); grid; xlabel('t/s'); ylabel('err/K'); title('Error'); legend('Error_{Hold}', 'Error_{Rate}');
figure;
subplot(3,1,1)
plot(t, heatRate_); grid; xlabel('t/s'); ylabel('Theta/°'); title('Aufheizrate');
plot(t, heatRate_); grid; xlabel('t/s'); ylabel('Theta/°'); title('Heat_{rate}');
subplot(3,1,2)
plot(t, P_rate_); grid; xlabel('t/s'); ylabel('Theta/°'); title('P_rate');
plot(t, y_rate_); grid; xlabel('t/s'); ylabel('Theta/°'); title('y_{rate}');
subplot(3,1,3)
plot(t, err_rate_); grid; xlabel('t/s'); ylabel('Theta/°'); title('Error');
plot(t, err_rate_); grid; xlabel('t/s'); ylabel('Theta/°'); title('Error_{rate}');
function [theta, s] = mass(si, dt, C, M, L, P, Td, dTheta)
s = si;
@@ -160,10 +226,10 @@ if ~isstruct(si)
if Td > 0
alpha = dt/Td;
end
s = struct('e', 0, 'a', alpha, 'x', 0, 'gain', 0.7);
s = struct('e', 0, 'a', alpha, 'x', 0, 'gain', 0.8);
end
s.e = s.e * (1-((L*dTheta)*dt)/(M*C));
s.e = s.e * (1-((L*(0+dTheta))*dt)/(M*C));
s.x = (1-s.a)*s.x + s.gain*s.a*P*dt;
s.e = s.e + s.x;
theta = s.e/(M*C);
@@ -194,6 +260,24 @@ end
endfunction
function [y, z] = isStable(zi, v, maxErr)
z = zi;
if ~isstruct(zi)
z = struct('v', 0, 'cnt', 3, 'cnt_reload', 3);
end
if z.cnt > 0
z.cnt = z.cnt - 1;
endif
if abs(v - z.v) <= maxErr
z.cnt = z.cnt_reload;
endif
y = (z.cnt == 0);
endfunction
function print_parameters(params)
align = 32;
names = fieldnames(params);