%Script for optimizing the drilling process with 4 inputs
%Needs external functions drilling4u.m, drillsolve4u.m 
%objfun4u.m and constraints.m.
%DEH - 2010

clear all
close all
clc

global db Wdbmax tt0 dstring Wdbt R0 K kc D dp Aa g rho rhos theta1

%Parameters:
db = 0.254; %Bit diameter, m
Wdbmax = 178.583; %Maximum WOB per m diameter, tons/m
tt0 = 10; %Trip time, hours
dstring = 0.1; %Diameter of drill string, m
D = 3000; %Depth, m
Aa = (db^2-dstring^2)*pi/4; %Cross-sectional area of annulus, m2
g = 9.81; %Gravity
R0 = 5; %Formation drillability, m/hr 
tempR0 = R0;
K = 15*R0; %Formation abrasiveness constant, hours 
kc = 5e5/R0; %Specific cutting force, N/m2-ton 
Wdbt = 63/R0; %Threshold WOB per m diameter, tons/m
theta1 = 900; %Friction parameter
rhos = 2700; %Formation density, kg/m3
rho = 1400; %Mud density, kg/m3
dp = 0.005; %Particle diameter, m

%Optimization:
u0 = [30;   %WOB, tons
    100;    %RPM
    40;     %qmud, l/s
    10];    %pc, bar

lb = [0; 0; 0; 0];
ub = [Wdbmax*db; 200; 50; 100];
options = optimset('Display','iter','TolFun',1e-12,'MaxFunEvals',300);
[u,fval,exitflag,output,lambda,grad,hessian] = fmincon(@objfun4u, u0,...
    [],[],[],[],lb,ub,@constraints,options)
Juu = hessian;

%Nominal optimum u and y:
y0 = drilling4u(u)