How to fit the data?

Below err is defined as the residual sum of squares. Using ?NumericQ on one of the arguments prevents FindMinimum from exact differentiation. In that way no problem arises even though the point with x == 0 is included.

err[a_, b_, c_, d_, e_?NumericQ] = Total[(a*Cosh[b*#^c*Sin[d*#^e]] - #2)^2 & @@@ data];

res[x_] = a*Cosh[b*x^c*Sin[d*x^e]] /. Last[FindMinimum[err[a, b, c, d, e], {
                            {a, 3, 1/100, 25}, {b, 2, 1/100, 25}, {c, 1, 1/100, 25},
                            {d, 23/10, 1/100, 25}, {e, 2, 1/100, 25}}, Method -> "InteriorPoint"]]

Show[Plot[res[x], {x, 0, 1.65}, PlotRange -> All], ListPlot[data]]

Alternative start values result in a smaller sum of squares:

FindMinimum[err[a, b, c, d, e],
  {{a, 2.6174, 1/100, 25}, {b, 1.7195, 1/100, 25}, {c, 2.3092, 1/100, 25},
   {d, 1.5033, 1/100, 25}, {e, 1.845, 1/100, 25}}, Method -> "InteriorPoint"]

{0.0011296543, {a -> 2.6174825, b -> 1.7194932, c -> 2.3092448, d -> 1.5033314, e -> 1.845972}}

These are highly likely optimal: When x == 0 the regression formula equals just a and because of the first data point a should presumably be close to 2.61.

If the last factor and Cosh are cancelled from the data values (using the right branch of ArcCosh) we get something we 2 peaks at high x-values:

data2 = Thread[{data[[All, 1]], PadLeft[{-1, -1}, 18, 1] ArcCosh[data[[All, 2]]/2.61]}]
ListPlot[data2]

We can solve for b and c such that the peaks are interpolated leaving only d and e for estimation:

sol = First[Solve[b #1^c Sin[d #1^e] == #2 & @@@ data2[[{-1, -6}]], {b, c}] /. C[1] -> 0 // Chop]

For different values of d and e let's look at the residual sum of squares excluding x == 0, because sol can't be evaluated in that instance

err0[a_, b_, c_, d_, e_?NumericQ] = Total[(a*Cosh[b*#^c*Sin[d*#^e]] - #2)^2 & @@@ Rest[data]];
search = Table[{d, e, Log[If[Im[#] == Im[#2] == 0, 
            Quiet[Min[err0[2.61, #, #2, d, e], 10^5.]], 10^5.]]} & @@
              ({b, c} /. sol), {d, 1/50, 7, 1/50}, {e, 1/50, 7, 1/50}] // Catenate;
ListPointPlot3D[search, PlotRange -> All, AxesLabel -> {x, y, z}]

the 2 smallest local minima of which coincide with the previous fits.

You can minimize the cost function by using NMinimize

  ClearAll["Global`*"]

data = {{0., 2.61}, {0.1, 2.62}, {0.2, 2.62}, {0.3, 2.62}, {0.4, 
2.63}, {0.5, 2.63}, {0.6, 2.74}, {0.7, 2.98}, {0.8, 3.66}, {0.9, 
5.04}, {1., 7.52}, {1.1, 10.74}, {1.2, 12.62}, {1.3, 10.17}, {1.4,
 5}, {1.5, 2.64}, {1.6, 11.5}, {1.65, 35.4}};

        model[x_] := a Cosh[b x^c Sin[d x^e]]
        cost = Total[(#2 - model[#1])^2 & @@@ data];
        fit = NMinimize[{cost, 1 < a < 3, 1 < b < 3, 1 < c < 3, 1 < d < 3, 1 < e < 3},
     {a, b, c, d, e}, Method -> "DifferentialEvolution"]

{0.00112965, {a -> 2.61748, b -> 1.71949, c -> 2.30924, d -> 1.50333, e -> 1.84597}}

Thread[{a, b, c, d, e} = {a, b, c, d, e} /. Last@fit];

Show[ListPlot[data], Plot[model[k], {k, 0, 1.7}]]

You can also use NonlinearModelFit with appropriate constrain on parameters and method choice.

ClearAll["Global`*"]

data = {{0., 2.61}, {0.1, 2.62}, {0.2, 2.62}, {0.3, 2.62}, {0.4, 
2.63}, {0.5, 2.63}, {0.6, 2.74}, {0.7, 2.98}, {0.8, 3.66}, {0.9, 
5.04}, {1., 7.52}, {1.1, 10.74}, {1.2, 12.62}, {1.3, 10.17}, {1.4,
 5}, {1.5, 2.64}, {1.6, 11.5}, {1.65, 35.4}};

model[x_] := a Cosh[b x^c Sin[d x^e]]
nlm = NonlinearModelFit[
  data, {model[x], 1 < a < 3, 1 < b < 3, 1 < c < 3, 1 < d < 3, 
   1 < e < 3}, {a, b, c, d, e}, x, 
  Method -> {"NMinimize", {Method -> "DifferentialEvolution"}}, MaxIterations -> 1000]

nlm // Normal

2.61749 Cosh[1.71949 x^2.30926 Sin[1.50334 x^1.84596]]

nlm["BestFitParameters"]

{a -> 2.61749, b -> 1.71949, c -> 2.30926, d -> 1.50334, e -> 1.84596}

Same picture

Just an extended comment that for this particular data/model combination it is difficult to find the values of parameters that minimize the residual sum of squares. I think there are two main issues:

1. The surface is extremely bumpy which can get any iterative search algorithms to get lost. (This is not necessarily a Mathematic, Maple, R, SAS, or *MathCAD" issue.)

2. Some of the parameter estimators are highly correlated with each other. That also causes problems for finding the optimal parameter values.

Below is a contour plot of the sum of squares evaluated at @user64494 's optimal values of $a$, $b$, and $c$ in the neighborhoods of $d$ and $e$. One can see the extreme bumpiness and the highly correlated nature of the estimators of $d$ and $e$. The red dot indicates the location of the optimal values of $d$ and $e$.