We optimize the objective function - with energy costs.

Input := 


Pcost[t_] = -.005 t -.5 asol[t] + 3.5 bsol[t] +

 .25 csol[t]/.solk[[2]]; 

Input := 


Plot[Pcost[t],{t,0,200}]

Output =


-Graphics-

Input := 


solPcost = FindRoot[Pcost'[t]==0,{t,40}]

Output =


{t -> 27.5171}

Input := 


Pcost[t]/.solPcost[[1]]

Output =


1.84894