function YWTwoCompExtendedLoading_FindOpt_simple()
%% Initial Simulation comparison (if on) 
    % Create fit object
	pf = ParameterFit();
    
% Pull in experimental data
	data = cell(1, 1);
    num = xlsread('200gL_8min_simple');
    data{1}=((num(1:end, 2))./1.51./.2/1000./153); %converting wavelength to concentration 
    times=num(1:end, 1).*60;                       %conver min to s

    % Create simulation
    sim = createSimulation(times, num);
    res=sim.run();

% Parameter Optimization 
% Set model parameters and enable sensitivities
    params = cell(6, 1);
    params{1}=makeSensitivity(0, {'MCL_KA'}, 0, -1, -1, 0, -1);
    params{2}=makeSensitivity(0, {'MCL_KD'}, 0, -1, -1, 0, -1);
    params{3}=makeSensitivity(0, {'MCL_QMAX'}, 0, -1, -1, 0, -1);
    params{4}=makeSensitivity(0, {'MCL_KA'}, 1, -1, -1, 0, -1);
    params{5}=makeSensitivity(0, {'MCL_KD'}, 1, -1, -1, 0, -1);
    params{6}=makeSensitivity(0, {'MCL_QMAX'}, 1, -1, -1, 0, -1);
	sim.setParameters(params, true(6, 1));

% Specify which components are observed in which factor for each
% observation / wavelength
	idxComp = cell(1, 1);
	idxComp{1} = [1, 1];

% Add the experiment to the fit (unit operation 0 is observed, which is the GRM) including name
% of the experiment and the different observations / wavelengths
	pf.addExperiment(data, sim, [0], idxComp, [], [], [], [], 'Composition Profiles', {'Total Protein'});

% Specify initial parameters and their lower and upper bounds
    initParams = [70, 1,3.6, 500, 0.01,5];  
	loBound = initParams.*0.1;
	upBound = initParams.*10;
    
% This variable serves as storage for the plot handles returned by the plot function of ParameterFit
% in the OutputFcn below.
	plotHd = [];

% Set some options (tolerances, enable Jacobian) and request some output (iteration statistics, plots)
	opts = optimset('TolFun', 1e-8, 'TolX', 1e-8, 'MaxIter', 100, 'Jacobian', 'on', 'Diagnostics', 'off', 'Display', 'iter', 'OutputFcn', @progressMonitor);
 
% Invoke optimizer lsqnonlin
	[optimalParams, optimalRes, ~, exitflag] = lsqnonlin(@residual, initParams, loBound, upBound, opts);
	success = (exitflag > 0);
    
     save('Set1-1.mat', 'optimalParams', 'optimalRes')

	function [varargout] = residual(x)
		%RESIDUAL Residual function that is passed to lsqnonlin, just forwards to ParameterFit object
		varargout = cell(nargout,1);
		[varargout{:}] = pf.residualVector(x);
	end

	function stop = progressMonitor(x, optimValues, state)
		%PROGRESSMONITOR Callback function for progress report invoked by lsqnonlin, just calls ParameterFit's plot function
		stop = false;
		if strcmp(state, 'iter')
			% Call plot function and reuse plot handles from previous call (storage outside this function in captured variable plotHd)
			plotHd = pf.plot(plotHd, optimValues.iteration, [optimValues.resnorm, optimValues.stepsize]);
		end
    end
end 


function sim = createSimulation(times, num)
	% General rate model
	mGrm = GeneralRateModel();

	mGrm.nComponents = 2; % Ordering: Protein, Arbitrary component 2 
	mGrm.nCellsColumn = 30; % Attention: This is very low and only used for illustration (short runtime)
	mGrm.nCellsParticle = 15; % Attention: This is very low and only used for illustration (short runtime)
	mGrm.nBoundStates = ones(mGrm.nComponents, 1);

	% Initial conditions
	mGrm.initialBulk = [0.0, 0.0];
	mGrm.initialSolid = [0.0, 0.0];
		
	% Transport
	mGrm.dispersionColumn          = 5.6e-8;
	mGrm.filmDiffusion             = [1e-6, 1e-6]; 
	mGrm.diffusionParticle         = [5e-11, 5e-11];
	mGrm.diffusionParticleSurface  = [0, 0];
    v = num(5,11);% Q = X mL/min;A = 0.34 cm2;epsilon = 0.35
	mGrm.interstitialVelocity      = v; 

	% Geometry
	mGrm.columnLength        = 0.059;
	mGrm.particleRadius      = 3.8e-5;
	mGrm.porosityColumn      = 0.34;
	mGrm.porosityParticle    = 0.91;

	% Adsorption
	mLangmuir=LangmuirBinding();
	mLangmuir.kineticBinding = true;
	mLangmuir.kA         = [7, 0.483];
	mLangmuir.kD         = [0.1, 0.001];
    mLangmuir.qMax       = [3.65, 4.2];
   	mGrm.bindingModel = mLangmuir;

	% Inlet
	mIn = PiecewiseCubicPolyInlet();
	mIn.nComponents = 2;
	
	mIn.constant       = zeros(2, mGrm.nComponents);
	mIn.linear         = zeros(2, mGrm.nComponents);
	mIn.quadratic      = zeros(2, mGrm.nComponents);
	mIn.cubic          = zeros(2, mGrm.nComponents);

	% Sec 1
    C = num(6,11);
    Perc_A =  num(7,11);
    Perc_B = num(8,11);
	mIn.constant(1,1)  = Perc_A*C/153;  
    mIn.constant(1, 2) = Perc_B*C/150;  
	
    % Sec 2
	% Everything is 0 already so no need to re-specify 
    
	% Model system
	mSys = ModelSystem();
	mSys.models = [mGrm, mIn];
	mSys.connectionStartSection = [0];
	mSys.connections = {[1, 0, -1, -1, -1, -1, 1.0]};

	% Configure simulator
	sim = Simulator.create();
	sim.sectionTimes = num(1:3, 3).*60;
	sim.sectionContinuity = false(1, 1);
	sim.solutionTimes = times;

	% Assign model
	sim.model = mSys;
end


% =============================================================================
%  CADET - The Chromatography Analysis and Design Toolkit
%  
%  Copyright (C) 2008-2020: The CADET Authors
%            Please see the AUTHORS and CONTRIBUTORS file.
%  
%  All rights reserved. This program and the accompanying materials
%  are made available under the terms of the GNU Public License v3.0 (or, at
%  your option, any later version) which accompanies this distribution, and
%  is available at http://www.gnu.org/licenses/gpl.html
% =============================================================================