%RSR_FIND_CONSTRAINTS

% Run fzero on rsr_optres to identify when different constraints switch
% between active and inactive for different values of disturbances

close all
clear all
%clc
load param.mat

% Constraint identifiers:
% 1: xB(2) > xB(2),min
% 2: Mr < Mr,max
% 3: Tr < Trmax
% 4: V < Vmax
% 5: F < Fmax
% 6: D > 0
% 7: B > 0
% 8: P > 0
% 9: RCY > 0
% 10: PD < PURGE_MAX

param.constraintID = 3 ; %Tell the file which constraint we want to check
param.options = optimset('Algorithm','interior-point','Display','off','TolCon',1e-5,'TolFun',1e-6,'TolX',1e-6,'MaxIter',1000,'MaxFunEvals',1e5) ; %Options for fmincon
%param.output = 'Opt_1' ; % To make the file reaseprec.m return both
% inequality and equality constraints, F < Fmax
param.output = 'Opt_2' ; % To make the file reaseprec.m return both inequality and equality constraints, F = Fmax 
param.Trmax = 390 ;
param.Mrmax = 11000 ;
param.Vmax = 30 ;
param.Fmax = 5 ;
param.pmax = 5 ; 
param.xfa = [1 0 0]' ;
param.prices = [1 0.025 -2 -0.5]' ;
param.variable = 'fmax' ; %Specify which variable is replaced by "d"
param.scalecon = [1e0 1e-3 1e-2 1e-1 1e-1]' ;

runwhat = 'optimize' ;
% runwhat = 'find' ;

switch runwhat
case('optimize')
    d = 5:0.5:6.5 ;
%     d =0.025:0.025:0.525;
    res = zeros(size(d));
    if length(d) == 1
        %close all
        param.options.Display = 'final' ;
        [res,un,J,con,l] = rsr_optres(d,param);
        %un = un.*param.scale_un ;
        make_output
    elseif length(d) > 1
        %close all
        cons = zeros(10,length(d));
        lambda = zeros(10,length(d));
        J = zeros(1,length(d));
        for i=1:length(d)
            disp(i)
            [res(i),un,J(i),con,l] = rsr_optres(d(i),param) ;
            cons(:,i) = con(1:10) ; % XB(2), Mrmax, Trmax, Vmax, Fmax, D>0, B>0, P>0, R > 0
            lambda(:,i) = l ;
        end
        for i=8
           figure(i),plot(d,cons(i,:),'b',d,lambda(i,:),'r');
%            figure(i),plot(d,cons(i,:),'b');
        end
        figure(11),plot(d,J);
        %figure(12+3),plot(d,res) ;
    end
case('find')
    %close all
    options = optimset('TolX',1e-5,'Display','iter') ; % Options for fsolve
    sol= fzero(@rsr_optres,[0.1 1],options,param);
end