The following code solves the finite deformation of the PneuNet. It consists of two parts:
Code 1: Solves for the strain parameters
Code 2: Reconstructs the 3D deformation curves
The Mooney-Rivlin constitutive model is used in this implementation.

Code 1: Solving the strain parameters

ClearAll["Global`*"]
Print["Mooney Rivlin"];
(*Structural and loading parameters*)
ϕ0=30 Pi/180;
L=66.8 10^-3;
W=10 10^-3;
n=15;
Wr= W/Sin[ϕ0];
t1=0.5 10^-3;
t2=0.75 10^-3;
t3=0.75 10^-3;
H1=1 10^-3;
H2=3 10^-3;
H3=9 10^-3;
H4=1 10^-3;
H5=1 10^-3;
P=50 10^3;

(*Strain assumptions*)
γ11=ϵ11+z k1;
γ22=ϵ22+z k2;
γ33=ϵ33;
ϵ33=1/((ϵ11+1)(ϵ22+1))-1;
λ1=γ11+1;
λ2=γ22+1;
λ3=γ33+1;

 (*Mooney-Rivlin Strain energy density function*)
I1=λ1^2+λ2^2+λ3^2;
I2=λ1^2 λ2^2+λ2^2 λ3^2+λ3^2 λ1^2;
c1=617111.114;
c2=-187159.625;
E3integrand=c1 (I1-3)+c2 (I2-3);

(* (* if use Neo-Hookean Strain energy density function, please use the following code*)
I1=λ1^2+λ2^2+λ3^2;
c10=193980.702;
E3integrand=c10 (I1-3);*)

(*Calculation of the elastic energy*)
E31=Simplify[(t1+t2+t3)Integrate[E3integrand,{z,0-H,H1-H},Assumptions->Im[z]==0]];E32=Simplify[(t1+t2)Integrate[E3integrand,{z,H1+H2-H,H1+H2+H5-H},Assumptions->Im[z]==0]];
E33=Simplify[t2 Integrate[E3integrand,{z,H1+H2+H5-H,H1+H2+H3-H},Assumptions->Im[z]==0]];
E34=Simplify[(t2+t3)Integrate[E3integrand,{z,H1+H2+H3-H,H1+H2+H3+H4-H},Assumptions->Im[z]==0]];
E3=Simplify[2n Wr(E31+E32+E33+E34)];

(*Calculation of the potential energy*)
f1=P H3 Wr;
f2=P H3 2 t3;
Hp=H1+H2+H3/2;
E1=n f1 2 t3 ((k1 Hp+ϵ11) ((Sin[ϕ-ϕ0])^2)+(k2 Hp+ϵ22) ((Cos[ϕ-ϕ0])^2));
E2=n f2((k1 Hp+ϵ11) ((Cos[ϕ-ϕ0])^2)+(k2 Hp+ϵ22) ((Sin[ϕ-ϕ0])^2))W/Sin[-ϕ0] Cos[-ϕ0];

(*Applying the minimum energy method*)
Δ=-E1-E2+E3;
exp1=Simplify[D[Δ,k1]];
exp2=Simplify[D[Δ,k2]];
exp3=Simplify[D[Δ,ϵ11]];
exp4=Simplify[D[Δ,ϵ22]];
exp5=ϵ22+H k2; (*the netural plane*)
exp6=Simplify[D[Δ,ϕ]];

(*Solve the equations using FindRoot. Alternatively,one may use NSolve,continuation methods,or other numerical techniques to solve it. Higher precision may also required. For computational efficiency,we assume
𝜙=ϕ0,as this approximation holds true for the current case.*)
ϕ=ϕ0;
sol=FindRoot[{exp1,exp2,exp3,exp4,exp5},{{k1,1 10^-3},{k2,1 10^-3},{ϵ11,1 10^-3},{ϵ22,1 10^-3},{H,1 10^-3}}]

(*Note the obtained results maynot be correct. Need lots of methods to verify the results is coorect. Here we just simply check the derivative of Δ with respect to all the parameters  is zero*)
{k1,k2,ϵ11,ϵ22,H}={k1,k2,ϵ11,ϵ22,H}/.sol
Print[{exp1, exp2, exp3, exp4, exp5, exp6}]

Code 2: Reconstruct the 3D curves

α=(Sqrt[k1^2 Cos[ϕ]^2+k2^2 Sin[ϕ]^2]);
β=k1 Cos[ϕ]^2+k2 Sin[ϕ]^2;
τ=(k1-k2)Sin[ϕ]Cos[ϕ];

X[s_]:=s ((Cos[ϕ])^2)-((k1 s β ((Cos[ϕ])^2))/(α^2))+((k1 β ((Cos[ϕ])^2) Sin[s α])/(α^3))+s ((Sin[ϕ])^2)-((k2 s β ((Sin[ϕ])^2))/(α^2))+((k2 β Sin[s α] ((Sin[ϕ])^2))/(α^3));
Y[s_]:=((τ (k1+k2+(k1-k2) Cos[2 ϕ]) ((s/2)-((Sin[s α])/(2 α))))/(α^2));
Z[s_]:=((Cos[s α] (k1+k2+(k1-k2) Cos[2 ϕ]))/(2 (α^2)))-((k1+k2+(k1-k2) Cos[2 ϕ])/(2 (α^2)));


r[t_]:={X[t],Y[t],Z[t]} ;
uT[t_]:=Simplify[r'[t]/Norm[r'[t]],t∈Reals]; (*the code here get the three vectors*)
vN[t_]:=Simplify[uT'[t]/Norm[uT'[t]],t∈Reals];
vB[t_]:=Simplify[Cross[r'[t],r''[t]]/Norm[Cross[r'[t],r''[t]]],t∈Reals];
vX={1,0,0} ;
vY={0,1,0};
vZ={0,0,1};
vT[s_]=-uT[s];
R2=1/k2;
b=R2 Sin[w/R2];
n=R2-R2 Cos[w/R2];
t=0;
Amatrix={{vX.vB[s], vX.vN[s], vX.vT[s]},{vY.vB[s], vY.vN[s], vY.vT[s]},{vZ.vB[s],vZ.vN[s],vZ.vT[s]}} ; (*Transformation matrix between two coordinate*)
Amatrix.{b,n,t};

shape=Simplify[{X[s],Y[s],Z[s]}+Amatrix.{b,n,t}]-h vN[s];  
shape1=Simplify[{X[s],Y[s],Z[s]}+Amatrix.{b,n,t}];

fig=ParametricPlot3D[{shape1[[1]],shape1[[2]],shape1[[3]]},{s,0,L},{w,-W/2,W/2},ColorFunction->Function[{x,y,z,t},Hue[t]],PlotStyle->Thick,AxesLabel->{"x","y","z"},BoxRatios->Automatic]