Numerical problems error

P_0 = 10e3;

kappa = 2 ;

Rfc = 3;

Nfc = 400;

Uoc = 5;

Pemin = -10e3;

Pemax = 75e3;

% Electric Motor

% Pi = d0 + d1*Pk + d2*Pk^2

d0 = 364.02;

d1 = 1.1014;

d2 = 1.028e-5;

Pkmin = -40e3;

Pkmax = 40e3;

% Battery

dEbmax = 0;

% Driving Cycle

m = 2e3;

Afcd = 1.2;

cr = 0.017;

rho = 1.25;

g = 9.81;

load('DrivingCycle')

N = length(v);

dt = 1;

theta = zeros(size(v));

Preq = (m*a.*v + m*g*sin(theta).*v + cr*m*g*cos(theta).*v...

+ rho/2*Afcd*v.^3)'; % Row vector

 

%% Convex Optmization: Optimal Control Problem

% Variables

deb = sdpvar(1,N+1);

pk = sdpvar(1,N);

pe = sdpvar(1,N);

pf = sdpvar(1,N);

pi = sdpvar(1,N);

t = sdpvar(1,N);

% if ICEOnOff

% b = binvar(1,N);

% else

b = 1;

% end

M = 100;

% Normalize

P0 = 1e4;

E0 = 2e5;

% Constraints

constr = [

(diff(deb)/dt == -pi*P0/E0) :'Battery Energy'

(pi >= (d0 + d1*pk*P0 + d2*(pk*P0).^2)/P0) :'Battery Power'

(pk >= Pkmin/P0) :'K-min'

(pk <= Pkmax/P0) :'K-max'

(pe >= Pemin/P0*b) :'E-min'

(pe <= Pemax/P0*b) :'E-max'

% (pf >= (e0 + e1*pe*P0 + e2*(pe*P0).^2)/P0 - M*(1-b)) :'Fuel'

% (pf >= pe + (P_0 + (kappa/2*Rfc)*((Nfc*(Uoc - kappa))+ ...

% sqrt((Nfc*(Uoc - kappa))^2 - 4*Nfc*Rfc*(pe + P_0))))) :'Fuel'

% (pf >= pe + (P_0 + (kappa/2*Rfc)*((Nfc*(Uoc - kappa)))+sqrt((Nfc*(Uoc - kappa))))) :'Fuel'

% (((pf - pe - P_0)*(2*Rfc/kappa)-(Nfc*(Uoc - kappa))).^2 - (Nfc*(Uoc - kappa))^2 - 4*Nfc*Rfc*(pe + P_0)>=0) :'Fuel'

(t>=0)

(t.*t - (Nfc*(Uoc - kappa))^2 - 4*Nfc*Rfc*(pe + P_0) <=0 )

(pf >= pe + P_0 + (kappa/2*Rfc)*((Nfc*(Uoc - kappa)) + t)) :'Fuel'

(pf >= 0) :'Fuel Positive'

(pe + pk == Preq/P0) :'Propulsion'

(deb(1) == 0) :'Initial'

(deb(end) >= -dEbmax/E0) :'Final'

];

% Objective

obj = sum(pf)*P0*dt/(N*P0);

params.BarHomogeneous = 1;

% Solver Options

options = sdpsettings('solver','gurobi');

% options.gurobi.QCPDual = 1;

% options.gurobi.TuneTimeLimit = 6000;