function Result = optimal_operation(x0, Agco, Avap, Aihx, Qvap, nair, u, d, sc);
Tgco_air = u(1);
Tvap_air = u(2);
disp('Find optimal operation')
aux.sc = sc;
aux.d = d;
aux.u = u;
f = 'shaft_work';
g = [];
H = [];
HessPattern = [];

xcom_L = [0; -inf; 0; 0]./sc.com.x;
xcom_U = [inf; inf; inf; inf]./sc.com.x;
xgco_L = [zeros(d.ngco,1); zeros(d.ngco,1); -inf*ones(d.ngco,1); zeros(d.ngco, 1); 0; 0; 0]./sc.gco.x;
xgco_U = [inf*ones(d.ngco,1); inf*ones(d.ngco,1); inf*ones(d.ngco,1); inf*ones(d.ngco, 1); inf; inf; inf]./sc.gco.x;
xval_L = [0; 0; -inf]./sc.val.x;
xval_U = [1; inf; inf]./sc.val.x;
xvap_L = [0; 0; 0; 0; 0; 0; 0]./sc.vap.x;
xvap_U = [inf; inf; inf; d.Vvap; inf; inf; inf]./sc.vap.x;
xihx_L = [zeros(d.nihx, 1); zeros(d.nihx, 1); zeros(d.nihx, 1); zeros(d.nihx, 1); zeros(d.nihx, 1); 0]./sc.ihx.x;
xihx_U = [inf*ones(d.nihx, 1); inf*ones(d.nihx, 1); inf*ones(d.nihx, 1); inf*ones(d.nihx, 1); inf*ones(d.nihx, 1); inf]./sc.ihx.x;
x_L = [xcom_L; xgco_L; xval_L; xvap_L; xihx_L];
x_U = [xcom_U; xgco_U; xval_U; xvap_U; xihx_U];
Name = 'CO2 optimal operation';
pSepFunc = [];
fLowBnd = [];
b_L = [
       Qvap;
       0;
       0;
       0;
       nair;
       ];
b_U = [
       Qvap;
       Agco; 
       Avap;
       Aihx;
       nair;
       ];
A = xlsread('A.xls','Operation');
SCALER = [];
for i=1:size(A,1)
  SCALER = [SCALER; sc.x'];
end 
A=A.*SCALER;
c = 'Heatpump';
dc = [];
d2c = [];
ConsPattern = [];
c_L = zeros(size(sc.res));
c_U = zeros(size(sc.res));
x_min = [];
x_max = [];
f_opt = [];
x_opt = [];

A*x0-b_L;
A*x0-b_U;
x0-x_L;
x0-x_U;

Prob = ProbDef;
Prob = conAssign(f, g, H, HessPattern, x_L, x_U, Name, x0, pSepFunc, fLowBnd, A, b_L, b_U, c, dc, d2c, ConsPattern, c_L, c_U, x_min, x_max, f_opt, x_opt); 
Prob.d=d;
Prob.sc=sc;
Prob.u = u;
Result = tomRun('minos', Prob);
disp(Result.ExitText)
x = Result.x_k;
A*x-b_L;
A*x-b_U;
xopt = x.*sc.x;
