Questions MBPC

[Leo Bakker]

Dear all,

I am new to CASADI and try to make a MBPC for indoor climate control in buildings. 

Please find my first attempt here:  https://github.com/leobakker/MBPC-CASADI/blob/master/MBPC_first_try.py

The MBPC controls the heating energy while minimizing thermal discomfort and energy use using a simple dynamic thermal model for building and heating system.  

I defined the following next steps:

next steps:
	-1- include vectors for boundary conditions (e.g. Tout, Tset, Qint Qsolar)
	-2- include multiple controlled actuators (e.g. heating + cooling + solar shading)
	-3- include temperature range (between min and max temperature) in penalty function
	-4- test multiple optimizers

I could not find any information helping me to realize the first two steps. 

Any pointers on all steps would be very welcomed and appreciated.

Thanks in advance, 

Leo Bakker 



[Joakim Børlum]

Hi Leo,

I would take a look at https://web.casadi.org/blog/ocp/ and try to redo your example using Opti (https://web.casadi.org/docs/#document-opti).

This way it is quite simple to introduce constraints on both states, inputs and outputs -- even as vectors.

Joakim

[Lukas Hoffmann]

Hi Leo,

I only had a glimpse at your code, but I assume by 'boundary conditions' you mean exogenous inputs. These can be passed to Casadi by defining additional symbolic variables. See the following minimal example in Matlab syntax:

T_in = SX.sym('T_in'); %system state indoor temperature
T_mass = SX.sym('T_mass'); %system state mass temperature

states = [T_in;T_mass];

u = SX.sym('u'); % ctrl

T_outdoor = SX.sym('T_outdoor'); % exogenous parameter 1 
T_set = SX.sym('T_set'); % exogenous parameter 2
params = [T_outdoor;T_set];

xdot1 = ((T_outdoor-T_in)/R01 + (T_mass - Tin)/R12 + Qin + u)/C1;
xdot2 = ((T_in-T_mass)/R01)/C2;

rhs = [xdot1;xdot2];

quad = (T_in-T_set)**2 + 0.001*u**2  ; %objective   

f = Function('f',{states,[u;params]},{rhs,quad}); %RK4 M = 4; % RK4 steps per interval                DT = t/N/M;                X0 = MX.sym('X0',2);                U = MX.sym('U',1);                PAR = MX.sym('PAR',2);                X = X0;                Q = 0;                for j=1:M                    [k1, k1_q] = f(X, [U;PAR]);                    [k2, k2_q] = f(X + DT/2 * k1, [U;PAR]);                    [k3, k3_q] = f(X + DT/2 * k2, [U;PAR]);                    [k4, k4_q] = f(X + DT * k3, [U;PAR]);                    X=X+DT/6*(k1 +2*k2 +2*k3 +k4);                    Q = Q + DT/6*(k1_q + 2*k2_q + 2*k3_q + k4_q);                 end                 obj.F = Function('F', {X0, [U;PAR]}, {X, Q}, {'x0','p'}, {'xf', 'qf'}); %etc etc  %% formulate NLP solver  nlp_prob = struct('f', J, 'x', w, 'g', g, 'p', p); casadi_solver = nlpsol('nlp_solver', 'ipopt', nlp_prob);             % then call the solver with actual parameter values sol = casadi_solver('x0',initialStates,'p',actualParameter,'lbx',lbw, 'ubx',ubw, 'lbg',lbg, 'ubg',ubg);