%TEST_ACCURACY_TEMP: File used to test how accurately we can solve MCHE
%equations and recycles using MATLAB solvers instead of Unisim's internal
%ones.

clear all
clc

global h
h=actxserver('UnisimDesign.Application');

load x0
load ua
load disturbances

%First, test how accurately we can solve for UA when specifying
%temperatures (using Unisim case MCHE_tempspecs.usc

lb_temp = [4240 ; 500 ; 101/3600 ; -97; -97 ; -97; -159.15 ; -159.15 ; 40/3600 ; 0/3600 ;  6/3600] ;
ub_temp = [4440 ; 537 ; 103/3600 ;  -82;  -82 ;  -82 ; -155.15 ; -155.15 ; 45/3600 ; 4/3600 ; 10/3600] ;
scale_x_temp = diag(2.^(fix(log2(ub_temp - lb_temp)))) ;

lb_temp = scale_x_temp\lb_temp ;
ub_temp = scale_x_temp\ub_temp ;
parameters.scale_x = scale_x_temp ;

x0_temp_scaled = scale_x_temp\x0_temp ;

x0 = [x0_temp_scaled(1:3) ; x0_temp_scaled(9:11)] ;
T0 = x0_temp_scaled(4:8) ;

parameters.solver = 'fsolve' ;
parameters.ua = ua ;

% Upper and lower bounds excluding temperature variables
lb = [lb_temp(1:3) ; lb_temp(9:11)] ;
ub = [ub_temp(1:3) ; ub_temp(9:11)] ;

% Use fsolve to find T for all x between lb and ub

tempfail = 0 ;
fsolve_options = optimset('Display','off','MaxIter',50,'TolFun',1e-8,'TolX',1e-8,'Algorithm','trust-region-reflective') ; % ,'Algorithm','levenberg-marquardt' 'trust-region-reflective'
flags =zeros(1,150);
for i=6
    disp(['Changing ' num2str(i) 'th element of x'])
    x = x0 ;
    xpert = linspace(lb(i),ub(i),25) ;
	for j=1:25
        x(i) = xpert(j) ;
        [T,fval,exitflag] = fsolve(@MCHE_UA_RES,T0,fsolve_options,x,d,parameters) ;
        ind = (i-1)*25 + j ;
        resnorm = fval'*fval ;
        flags(ind) = exitflag ;
        if resnorm < fsolve_options.TolFun
            disp('Solution found')
        elseif resnorm > fsolve_options.TolFun
            if exitflag == 0
                disp('Maximum iterations used')
            elseif exitflag ~= 0
                disp('Did not converge')
            end
            tempfail = tempfail + 1 ;
        end
	end
end



% Use fmincon, with a dummy objective function.

% parameters.solver = 'fmincon' ;
% fmincon_options = optimset('Algorithm','interior-point','TolCon',1e-4,'Display','off','MaxIter',50,'MaxFunEvals',1000); % % 'interior-point' 'active-set'
% tempfail = 0 ;
% for i=6
%     disp(['Changing ' num2str(i) 'th element of x'])
%     x = x0 ;
%     xpert = linspace(lb(i),ub(i),25) ;
%     for j=1:25
%         x(i) = xpert(j) ;
% %         [T,fval,exitflag] = fmincon(@dummyJ,T0,[],[],[],[],[],[],@MCHE_UA_RES,fmincon_options,x,d,parameters) ;
%         [T,fval,exitflag] = fmincon(@dummyJ,T0,[],[],[],[],lb_temp(4:8),ub_temp(4:8),@MCHE_UA_RES,fmincon_options,x,d,parameters) ;
%         if exitflag > 0
%             disp('Solution found')
%         elseif exitflag == 0
%             disp('Maximum iterations used')
%         elseif exitflag < 0
%             disp('Did not converge')
%             tempfail = tempfail + 1 ;
%         end
%     end
% end

%Repeat, using lsqnonlin.

% lsqnonlin_options = optimset('Display','final','MaxIter',100,'TolFun',1e-8,'Algorithm','levenberg-marquardt') ; %
% parameters.solver = 'lsqnonlin' ;
% tempfail = 0 ;
% flags =zeros(1,150);
% for i=1:6
%     disp(['Changing ' num2str(i) 'th element of x'])
%     x = x0;
%     xpert = linspace(lb(i),ub(i),25) ;
%     for j=1:25
%         x(i) = xpert(j) ;
%         [T,resnorm,residual,exitflag] = lsqnonlin(@MCHE_UA_RES,T0,[],[],lsqnonlin_options,x,d,parameters) ;
%         disp(resnorm)
%         ind = (i-1)*25 + j ;
%         flags(ind) = exitflag ;
%         if resnorm < lsqnonlin_options.TolFun
%             disp('Solution found')
%         elseif resnorm > lsqnonlin_options.TolFun
%             if exitflag == 0
%                 disp('Maximum iterations used')
%             elseif exitflag ~= 0
%                 disp('Did not converge')
%             end
%             tempfail = tempfail + 1 ;
%         end
%     end
% end

%Reset the Unisim model to nominal values
Reset_MCHE_temp(x0_temp,d,h) ;