restart;

with(Statistics):with(plots):with(Optimization):with(LinearAlgebra):

# given data from strain rate curve
E_0[theta] := 7.883352314*10^9;
alpha[theta]:= 0.982

# experimental creep data under 44 at 100 degree celcius
c_strain := Vector ([<<0>,<0.0284698>,<0.0533808>,<0.0782918>,<0.0996441>,<0.124555>,<0.142349>,<0.156584>,<0.16726>,<0.177936>,<0.181495>,<0.188612>,<0.192171>,<0.19573>,<0.19573>,<0.202847>,<0.206406>,<0.206406>,<0.209964>,<0.209964>,<0.209964>,<0.206406>,<0.209964>>]):

c_time := Vector ([<<0>,<0>,<0.048>,<0.192>,<0.352>,<0.544>,<0.704>,<0.896>,<1.088>,<1.312>,<1.52>,<1.76>,<1.984>,<2.208>,<2.464>,<2.736>,<3.088>,<3.392>,<3.664>,<4.016>,<4.352>,<4.592>,<4.832>>]):
sigma[0] := 44*10^6;
epsilon[0] := sigma[0]/E_0[theta];

# change vector to list
c_strain := convert(c_strain,list):
c_time := convert(c_time,list):

# extract zero from list
c_strain := c_strain [2..-1];
c_time := c_time [2..-1];

# for further calculation need to know how many elements are in the list
M := nops(c_strain);
N := nops(c_time);

# constitutive equation
creep_strain := proc(t)
local i;
options operator, arrow;
epsilon[0]*(1 + alpha[theta]*add(-(B[i]*(-beta[i]*t + exp(-beta[i] *t) -1))/beta[i],i=1..3))
end proc;

#define objective function

Digits := 9:
obj := add(
         (
           creep_strain(c_time[j])
           -
           c_strain[j]
         )^2
         , j=1..N
       ):

# curve fitting
opt := NLPSolve(obj,
{
  beta[1] >= 0,
    beta[2] >= 0,
    beta[3] >= 0,
    B[1] >= 0, B[1]<= 0.2,
    B[2] >= 0, B[2] <= 0.3,
    B[3] >= 0, B[3] <= 0.4
}
);

# plot the result
display(
   ScatterPlot(c_time, c_strain, symbol=circle, color=blue),
 
  plot(eval(creep_strain(t), opt[2]), t=0.00..max(c_time), color=black)
)