Contents
Energy-Aware MPC Demo
This demo illustrates the usage of the EAMPC class on a 2-states, 1-input linear system subject to state and output noise. Several sensor nodes are used for disturbance rejection purposes. We assume that every node measures the full state vector, corrupted by an additive disturbance. At every time step, once every node got a measurement, an estimated state value is obtained by taking the mean value of the outputs. This estimation is transmitted to the controller in accordance with a threshold-based policy, where the estimated state is compared with a predicted value which has been precomputed by the controller and transmitted beforehand to the sensors. This policy is inteded to reduce the number of communications between controller and sensors, hence saving sensor nodes battery. In this example the energy-aware MPC is compared with a traditional MPC scheme, where the measurements are simply transmitted to the controller at every time step, and no predictions are computed and transmitted by the controller.
by D. Bernardini, 2010.
clear all close all
Plant model
State-space matrices of the discrete-time LTI system
Plant.A = [.21 -.39; -.39 .82]; Plant.B = [0 1]';
Network model
specify parameters for sensor nodes
Net.nodes = 3; % number of wireless sensor nodes Net.th = .05*ones(2,1); % transmission thresholds Net.Ne = 10; % estimation horizon
Controller design
% element-wise constraints on output, input, and input rate Limits.ymin = [-2 -2]'; % min output Limits.ymax = [2 2]'; % max output Limits.umin = -1; % min input Limits.umax = 1; % max input Limits.dumin = -.6; % min input rate Limits.dumax = .6; % max input rate % weights matrices for MPC objective function Weights.Qu = .1; % input weight Weights.Qy = eye(size(Plant.A,1)); % output weight Weights.Qn = Weights.Qy; % terminal weight Weights.rho = Inf; % positive weight for soft output constraints % (if rho=Inf then hard constraints are imposed) % other parameters Params.pnorm = 2; % norm used in the objective function (1, 2 or Inf). Params.N = 10; % prediction horizon Params.Nc = Params.N; % control horizon % build EAMPC object EAobj = eampc(Plant,Net,Limits,Weights,Params); % plot explicit MPC partition Options = mpt_options; Options.shade = .5; Options.samecolors = 1; figure plot(EAobj.Ctrl.Pn,Options); xlabel 'x_1' ylabel 'x_2' zlabel 'u_{prev}' axis(1.4*[Limits.ymin(1) Limits.ymax(1) ... Limits.ymin(2) Limits.ymax(2) Limits.umin Limits.umax]); title('Explicit MPC partition') box on
looking for available solvers... MPT toolbox 2.6.2 initialized... Copyright (C) 2003-2006 by M. Kvasnica, P. Grieder and M. Baotic Send bug reports, questions or comments to mpt@control.ee.ethz.ch For news, visit the MPT web page at http://control.ee.ethz.ch/~mpt/ LP solver: CDD Criss-Cross QP solver: quadprog MILP solver: CPLEX9 MIQP solver: CPLEX9 Vertex enumeration: CDD Run 'mpt_studio' to start the GUI. Run 'mpt_setup' to set global parameters. Region: 20/31 Region: 40/44 mpt_mpqp: 45 regions -> Generated 1 partition. Solution consists of 45 regions 45 regions with 7 different control laws control law # 1, 1 region --> 1 region control law # 2, 1 region --> 1 region control law # 3, 1 region --> 1 region control law # 4, 16 regions --> 7 regions control law # 5, 16 regions --> 7 regions control law # 6, 5 regions --> 1 region control law # 7, 5 regions --> 1 region controller partition reduced to 19 regions Plotting...

Simulation: energy-aware MPC vs standard MPC
The energy-aware MPC controller is compared with a standard MPC control scheme, where at every time step measurements are transmitted to the controller regardless to the threshold logic, and no predictions are transmitted to the sensors.
The two controllers are simulated in closed-loop with the LTI system
where Moreover, measurements are affected by an additive error, i.e.,
for all Net.nodes, where
is the output measured by the i-th node, and
, for all
.
% build EAMPC object for standard MPC MPCobj = EAobj; % copy EAMPC object MPCobj.Net.th = zeros(MPCobj.ny,1); % set thresholds to zero % simulation length Tsteps = 20; % disturbance realizations loadNoise = 1; % 0: pick random noise, 1: load noise data form Noise.mat if loadNoise==0 % additive disturbance bounds wmax = [.03 .03]'; % output noise bounds for each output and each node vmax = kron([.06 .06]',ones(1,Net.nodes)); % additive disturbance realizations W = diag(2*wmax)*rand(EAobj.ny,Tsteps)-wmax*ones(1,Tsteps); % output noise realizations V = zeros(EAobj.ny,EAobj.Net.nodes,Tsteps); for j=1:EAobj.Net.nodes V(:,j,:) = diag(2*vmax(:,j))*rand(EAobj.ny,Tsteps) ... -kron(ones(1,Tsteps),vmax(:,j)); end % initial state (randomly picked in the feasible state set) vert = extreme(projection(EAobj.Ctrl.Pfinal,1:EAobj.ny)); comb = rand(size(vert,1),1); comb = comb/sum(comb); x0 = (comb'*vert)'; save Noise wmax vmax W V x0 else load Noise.mat end % init energy-aware MPC simluation variables X = x0; % state Xestim = []; % estimated state U = []; % input Uprev = zeros(EAobj.nu,1); % previous input J = []; % experimental cost function EAobj = EAobj.init_sim(); % init simulation data % init standard MPC simluation variables XX = x0; % state XXestim = []; % estimated state UU = []; % input UUprev = zeros(MPCobj.nu,1); % previous input JJ = []; % experimental cost function MPCobj = MPCobj.init_sim(); % init simulation data % compute first sequence of output predictions and transmit them to sensors % (we assume that the initial state is known) EAobj = EAobj.send_predictions(X,Uprev); MPCobj = MPCobj.send_predictions(XX,UUprev); % note: since MPCobj.Net.th = 0, predictions are not actually transmitted here % run simulation for k=1,2,...,Tsteps for k=1:Tsteps % measurements acquisition [xestim,EAobj] = EAobj.get_measurements(X(:,end),V(:,:,k)); Xestim(:,end+1) = xestim; [xestim,MPCobj] = MPCobj.get_measurements(XX(:,end),V(:,:,k)); XXestim(:,end+1) = xestim; % compute the control law U(:,end+1) = EAobj.get_input(Xestim(:,end),Uprev); UU(:,end+1) = MPCobj.get_input(XXestim(:,end),UUprev); % if needed, compute new predictions EAobj = EAobj.send_predictions(Xestim(:,end),Uprev); MPCobj = MPCobj.send_predictions(XXestim(:,end),UUprev); % evaluate evperimental costs if EAobj.Params.pnorm == 2 J(k) = X(:,end)'*Weights.Qy*X(:,end) + ... U(:,end)'*Weights.Qu*U(:,end); else J(k) = norm(Weights.Qy*X(:,end),EAobj.Params.pnorm) + ... norm(Weights.Qu*U(:,end),EAobj.Params.pnorm); end if MPCobj.Params.pnorm == 2 JJ(k) = XX(:,end)'*Weights.Qy*XX(:,end) + ... UU(:,end)'*Weights.Qu*UU(:,end); else JJ(k) = norm(Weights.Qy*XX(:,end),MPCobj.Params.pnorm) + ... norm(Weights.Qu*UU(:,end),MPCobj.Params.pnorm); end % state evolution X(:,end+1) = Plant.A*X(:,end) + Plant.B*U(:,end) + W(:,k); XX(:,end+1) = Plant.A*XX(:,end) + Plant.B*UU(:,end) + W(:,k); % update previous input Uprev = U(:,end); UUprev = UU(:,end); end
Results
% plot state trajectory figure; subplot(2,1,1); plot(1:Tsteps,X(1,1:Tsteps),'b',1:Tsteps,XX(1,1:Tsteps),'r--','LineWidth',2); title('Energy-aware vs. standard MPC: state 1') ylabel 'x1' legend('EA-MPC','std. MPC','Location','Best') grid on subplot(2,1,2); plot(1:Tsteps,X(2,1:Tsteps),'b',1:Tsteps,XX(2,1:Tsteps),'r--','LineWidth',2); title('Energy-aware vs. standard MPC: state 2') ylabel 'x2' legend('EA-MPC','std. MPC','Location','Best') grid on % plot input trajectory figure stairs(1:Tsteps,U(1,1:Tsteps),'b','LineWidth',2) hold on stairs(1:Tsteps,UU(1,1:Tsteps),'r--','LineWidth',2) title('Energy-aware vs. standard MPC: control move') ylabel 'u1' legend('EA-MPC','std. MPC','Location','Best') grid on % plot WSN activity figure stairs(0:Tsteps,[0 EAobj.Sim.tx(1:Tsteps)]+1.5,'b','LineWidth',1.5); hold on stairs(0:Tsteps,[0 EAobj.Sim.rx(1:Tsteps)],'r','LineWidth',1.5); set(gca,'Ytick',[0 1 1.5 2.5],'YTickLabel',{'0' '1' '0' '1'}); title('Sensor-to-controller transmissions: time plots') xlabel 'Sampling instants' ylabel 'Transmissions' legend('out','in','Location','Best') grid on % plot pie graph for transmission data figure tx_perc = 100*sum(EAobj.Sim.tx)/Tsteps; rx_perc = 100*sum(EAobj.Sim.rx)/Tsteps; free_perc = 100-tx_perc-rx_perc; pie([tx_perc free_perc rx_perc],[0 1 0], ... {['Outgoing: ' num2str(tx_perc,'%3.1f') '%'], ... ['No transmissions: ' num2str(free_perc,'%3.1f') '%'], ... ['Incoming: ' num2str(rx_perc,'%3.1f') '%']}); title('Sensor-to-controller transmissions: overall statistics') % print experimental costs sumJ = sum(J); sumJJ = sum(JJ); fprintf('\n== Results for Energy-Aware MPC:\n'); fprintf('threshold TH = [%.2f %.2f]'', experimental cost J = %.3f.\n', ... EAobj.Net.th',sumJ); fprintf('\n== Results for standard MPC:\n'); fprintf('threshold TH = [%.2f %.2f]'', experimental cost J = %.3f.\n', ... MPCobj.Net.th',sumJJ); fprintf('\n== Energy-Aware vs standard MPC: transmission savings = %.2f%%\n', ... free_perc); fprintf('\n== Energy-Aware vs standard MPC: performance ratio = %.2f%%\n', ... (sumJ/sumJJ-1)*100);
== Results for Energy-Aware MPC: threshold TH = [0.05 0.05]', experimental cost J = 0.441. == Results for standard MPC: threshold TH = [0.00 0.00]', experimental cost J = 0.432. == Energy-Aware vs standard MPC: transmission savings = 60.00% == Energy-Aware vs standard MPC: performance ratio = 1.98%




Conclusions
In the proposed example, the energy-aware MPC scheme provided a reduction of 60% in the number of transmissions between the controller and the sensor nodes. This energy savings have been obtained with a small loss in the performance index (around 2%), with respect to a traditional MPC scheme. Indeed, the threshold values Net.th need to be properly tuned, especially in function of the magnitude of the disturbances, in order to obtain a valuable reduction in the transmission rate and also achieve a satisfying performance.