%% Cost_function.m
% Summary:      Insert and exstract data from unisim, can also with small 
%               modifications be used as the cost functin in an 
%               optimizing algorithm
% author:       Tor Anders Marvik 
% organization: Department of Chemical Engineering, NTNU, Norway
% contact:      torandma@stud.ntnu.no  
% requires:     matlab and unisim
% todo 1.0:     All the data should advantageously be orgenized in a 
%               structure array. 
% started       Oktober 2012
%Input:         setpoints to unisim 
%Output:        vectores and stucture array with unisim data
%%

function [cost_parameters,merke,c,Temp_profile,data_all] = cost_function(x,hyCase,data_all,i)
stream=hyCase.Flowsheet.MaterialStreams;                    %stream path
EnergyStreams=hyCase.Flowsheet.EnergyStreams;               % energy strems 
operations=hyCase.Flowsheet.Operations;
col=hyCase.Flowsheet.Operations.Item('21_VE103');
sol=hyCase.Solver;
script_folder='C:\Users\torandma\Documents\prosjekt karsto\matlab';
merke=0;
in=0;

%Setting controller setpoints and other degrees of freedom
sol.canSolve=0;
stream.Item('1b').TemperatureValue=x(1);  % C   | Temp before feed sep
stream.Item('13').PressureValue=x(2);         % kPa | Pressure feed sep 
operations.Item('21_PV2001').PressureDrop.Value=x(3); %kPa pressure drop cold sep. - column
RC(1-x(4),script_folder,hyCase)                  % frac| Split ratio liq. cold sep  
EnergyStreams.Item('Q_side_btm').HeatFlowValue=x(5)/3600; %kj/h | energy bmt side boiler
col.ColumnFlowsheet.Specifications.Item('Temperature').GoalValue=x(6); %setpoint column bottom temperature 
stream.Item('14').PressureValue=x(7);    %kPa |Pressure after expander
stream.Item('D').PressureValue=x(7)-200;
stream.Item('12').PressureValue=x(7)-200; % kPa | Pressure diff cold sep
sol.canSolve=1;
  
isknown=0;   % parameters to check if unisim converges
isknown=isknown + stream.Item('D').MassFlow.isKnown         -1;
isknown=isknown + stream.Item('B').MassFlow.isKnown         -1;
isknown=isknown + EnergyStreams.Item('E1').HeatFlow.isKnown -1;
isknown=isknown + stream.Item('1').Temperature.IsKnown      -1;
isknown=isknown + stream.Item('1').Temperature.IsKnown      -1;
isknown=isknown + stream.Item('1').Temperature.IsKnown      -1;
isknown=isknown + stream.Item('1b').Temperature.IsKnown     -1;

isknown=isknown + strcmp('Converged',operations.Item('RCY-7').RecycleConvergence) -1;
isknown=isknown + strcmp('Converged',operations.Item('RCY-5').RecycleConvergence) -1;
isknown=isknown + strcmp('Converged',operations.Item('RCY-1').RecycleConvergence) -1;
isknown=isknown + strcmp('Converged',operations.Item('RCY-6').RecycleConvergence) -1;

if isknown<0
    fprintf('Unisim failed to converge\n')
    fprintf('Please solve the problem manualy\n\n')
    in = input('Has the problem been solved manualy? (y=0/n=1):');   
end


if in == 1                      %model not converged, insert zero vector
    cost_parameters=zeros(1,5);
    Temp_profile=zeros(1,39);
    c=zeros(1,8);
    merke=1;
else                            %model converged, extract data
    

    %===========cost function parameters================================
    D=stream.Item('D').MassFlowValue*3600;  %kg/h
    B=stream.Item('B').MassFlowValue*3600;    % kg/h
    V=col.ColumnFlowsheet.MaterialStreams.Item('Boilup').MassFlowValue*3600;    %reboiler energy

    powercomp0=operations.Item('K-100').energyValue*3600;  %kJ/h | duty sales gass comp
    powercomp1=operations.Item('K-101').energyValue*3600;  %kJ/h | duty sales gass comp
    powercomp2=operations.Item('K-102').energyValue*3600;  %kJ/h | duty sales gass comp
    powercomp3=operations.Item('K-103').energyValue*3600;  %kJ/h | duty sales gass comp
    powercomp4=operations.Item('K-104').energyValue*3600;  %kJ/h | duty sales gass comp
    powercomp5=operations.Item('K-105').energyValue*3600;  %kJ/h | duty sales gass comp

    comp_salesgas=powercomp1+powercomp1+powercomp2;
    comp_boster=powercomp3+powercomp4+powercomp5;
    %==================================================================

    %=================constraints======================================
    c(1,1)=col.ColumnFlowsheet.PressuresValue(end)/100 ;      %[bar] Tower btm pressure
    c(1,2)=stream.Item('13').PressureValue/100;                %[bar] expansion pressure
    c(1,3)=col.columnFlowsheet.temperaturesValue(20) ;    %Temperature stage 19
    c(1,4)=operations.Item('21_PV2001').PressureDropValue/100; %[bar] pressure drop column cold sep
    c(1,5)=stream.Item('20').TemperatureValue;                %[C] Temperature befoure compressor
    c(1,6)=stream.Item('B').componentMolarFractionValue(4);     %methane bottom
    c(1,7)=stream.Item('D').componentMolarFractionValue(3);     %CO2 top
    c(1,8)=operations.Item('21_HA204').UAValue*3600;       %UA kJ/C*h

    Temp_profile=col.columnFlowsheet.TemperaturesValue;

    merke=0;
    cost_parameters=[D, B, V, comp_salesgas, comp_boster];
    %cost=[-1,-2,1,1,1]*cost_parameter,s';   %coment out if the function is to be used
                                             %in the automized optimization

    %Creating a struct of all important variables
    data_all.constraints.tower_btm_p(:,i)=col.ColumnFlowsheet.PressuresValue(end)/100;
    data_all.constraints.expansion_pressure(:,i)=stream.Item('13').PressureValue/100;
    data_all.constraints.temp_tray19(:,i)=col.columnFlowsheet.temperaturesValue(20);
    data_all.constraints.dela_p_column_coldsep(:,i)=operations.Item('21_PV2001').PressureDropValue/100;
    data_all.constraints.temp_before_comp(:,i)=stream.Item('20').TemperatureValue;
    data_all.constraints.C1_btm(:,i)=stream.Item('B').componentMolarFractionValue(4);
    data_all.constraints.CO2_top(:,i)=stream.Item('D').componentMolarFractionValue(3);
    data_all.constraints.UA_21HA204(:,i)=operations.Item('21_HA204').UAValue*3600;

    data_all.cost_param.D_flow(:,i)=stream.Item('D').MassFlowValue*3600;
    data_all.cost_param.B_flow(:,i)=stream.Item('B').MassFlowValue*3600;    % kg/h
    data_all.cost_param.V_flow(:,i)=col.ColumnFlowsheet.MaterialStreams.Item('Boilup').MassFlowValue*3600;
    data_all.cost_param.compwork_booster(:,i)=comp_boster;
    data_all.cost_param.compwork_salesgas(:,i)=comp_salesgas;

    %stream data
    data_all.feed_sec2=extract_stream_data('feed_sec2',i,stream,data_all.feed_sec2);
    data_all.s1=extract_stream_data('1',i,stream,data_all.s1);
    data_all.s1b=extract_stream_data('1b',i,stream,data_all.s1b);
    data_all.s2=extract_stream_data('2',i,stream,data_all.s2);
    data_all.s3=extract_stream_data('3',i,stream,data_all.s3);
    data_all.s4=extract_stream_data('4',i,stream,data_all.s4);
    data_all.s5=extract_stream_data('5',i,stream,data_all.s5);
    data_all.s6=extract_stream_data('6',i,stream,data_all.s6);
    data_all.s7=extract_stream_data('7',i,stream,data_all.s7);
    data_all.s8=extract_stream_data('8',i,stream,data_all.s8);
    data_all.s9=extract_stream_data('9',i,stream,data_all.s9);
    data_all.s10=extract_stream_data('10',i,stream,data_all.s10);
    data_all.s11=extract_stream_data('11',i,stream,data_all.s11);
    data_all.s12=extract_stream_data('12',i,stream,data_all.s12);
    data_all.s13=extract_stream_data('13',i,stream,data_all.s13);
    data_all.s14=extract_stream_data('14',i,stream,data_all.s14);
    data_all.s15=extract_stream_data('15',i,stream,data_all.s15);
    data_all.s16=extract_stream_data('16',i,stream,data_all.s16);
    data_all.s17=extract_stream_data('17',i,stream,data_all.s17);
    data_all.s18=extract_stream_data('18',i,stream,data_all.s18);
    data_all.s19=extract_stream_data('19',i,stream,data_all.s19);
    data_all.to_sec7=extract_stream_data('to_sec7',i,stream,data_all.to_sec7);
    data_all.D=extract_stream_data('D',i,stream,data_all.D);
    data_all.B=extract_stream_data('B',i,stream,data_all.B);
    data_all.side_btm_out=extract_stream_data('side_btm_out',i,stream,data_all.side_btm_out);
    data_all.side_btm_in=extract_stream_data('side_btm_in',i,stream,data_all.side_btm_in);
    data_all.side_top_out=extract_stream_data('side_top_out',i,stream,data_all.side_top_out);
    data_all.side_top_in=extract_stream_data('side_top_in',i,stream,data_all.side_top_in);


    

end

end













