The 4-peg tower is significantly more difficult than the original version with 3 pegs. I use Jaime Rangel-Mondragon’s code for computation and add simple rules for geometric transformation to create the animation for moves.

Notation

Run the code in the “Jaime Rangel-Mondragon’s code” section below to load the definition for the function HanoiMoves. In my example, I put 10 plates in the game with 4 pegs or needles.
In[]:=
allConfigs=HanoiMoves[10,4];
Say we extract the first 5 steps. The columns in the matrix stands for needle or peg:
In[]:=
allConfigs[[;;5]]//MatrixForm
Out[]//MatrixForm=
{1,2,3,4,5,6,7,8,9,10}
{}
{}
{}
{2,3,4,5,6,7,8,9,10}
{1}
{}
{}
{3,4,5,6,7,8,9,10}
{1}
{2}
{}
{4,5,6,7,8,9,10}
{1}
{2}
{3}
{4,5,6,7,8,9,10}
{1}
{}
{2,3}
The list of number in each column stands for the pile of numbered plated sitting at the current location. For instance the last row in the above matrix has
{{4,5,6,7,8,9,10},{1},{},{2,3}}
, I can illustrate it as the following diagram:

Use GeometricTransformation for animation

Define two association to record the initial location and position of the center for each plate at the left most pile with
n
plates:
In[]:=
plates[n_Integer]:=AssociationThread[Range[n]->Array[Style[Rectangle[{-.45#/n,1/n*(n-#)},{.45#/n,1/n*(n-#+1)},RoundingRadius->.3/n],Hue[#/n+0.2]]&,n]​​]
In[]:=
platesCenterInit[n_]:=AssociationThread[Range[n]->Array[{0,1./n*(n-#)+1/(2n)}&,n]]
For 10 plates:
In[]:=
data1=plates[10];​​data2=platesCenterInit[10];​​Graphics[{Line[{{0,-0.05},{0,1.1}}],Values@data1,Point[Values[data2]]},AspectRatio->1.2]
Out[]=
Given any entry of the move set from Jaime’s function, the new target location for each plate is computed with following rule:
pos[plts_List,needleId_List,n_Integer]:=Module[{len=Length[plts],id=needleId[[1]]},​​If[plts!={},​​AssociationThread[plts->Array[{id-1,(len-#)/n+1/(2n)}&,len]],​​<||>​​]]
targetAssoc[config_List,n_Integer]:=Join@@MapIndexed[pos[#1,#2,n]&,config]
For instance at Step=35
In[]:=
allConfigs[[35]]
Out[]=
{{},{3,4,5,6},{2},{1,7,8,9,10}}
The new location for each center of the plates become:
In[]:=
target=targetAssoc[%,10]
Out[]=
31,
7
20
,41,
1
4
,51,
3
20
,61,
1
20
,22,
1
20
,13,
9
20
,73,
7
20
,83,
1
4
,93,
3
20
,103,
1
20

Move plates to new locations with GeometricTransformation:
In[]:=
gConfig2=Table[GeometricTransformation[data1[i],TranslationTransform[target[i]-data2[i]]],{i,10}];
In[]:=
move=allConfigs[[35]];​​pegs=Array[Line[{{#,-0.2},{#,1.1}}]&,4,0];​​ground=Line[{{-.45,0},{3.45,0}}];​​Graphics[Table[Inset[Style[Length[move[[i]]],14],{i-1,1.25}],{i,4}]~Join~{​​pegs,ground,gConfig2​​},​​AspectRatio->.5,PlotRange->{{-0.45,3.45},{-.05,1.3}}]
Out[]=
I compile all animation code in one place in the “Utility- Animation code” section below. Just run that part you should be fine to create the animation for all moves.

OEIS A007664

You can compute the steps directly with this WL function used on the OEIS page:
In[]:=
a[n_]:=a[n]=Min[Table[2*a[k]+2^(n-k)-1,{k,0,n-1}]];a[0]=0;Table[a[n],{n,0,48}](*Jean-FrançoisAlcover,Dec06×2011,afterM.F.Hasler*)
Out[]=
{0,1,3,5,9,13,17,25,33,41,49,65,81,97,113,129,161,193,225,257,289,321,385,449,513,577,641,705,769,897,1025,1153,1281,1409,1537,1665,1793,2049,2305,2561,2817,3073,3329,3585,3841,4097,4609,5121,5633}

Utility- Jaime Rangel-Mondragon’s code

In[]:=
sums[s_,i_]:={}/;(s<i)||(i==0)​​sums[s_,s_]:={Table[1,{s}]}​​sums[s_,1]:={{s}}​​sums[s_,i_]:=Module[{d},​​Flatten[Table[Map[Join[{d},#]&,sums[s-d,i-1]],{d,s-i+1}],1]]
In[]:=
hanoiP[n_,p_]:=Join[Table[1,{n-1}],Table[0,{p-n-1}],{2n-1}]/;n<p-1​​hanoiP[n_,p_]:=hanoiP[n,p]=Module[{v,t},​​t=sums[n-1,p-2];​​v=Map[Join[#,{2Plus@@Table[Last[hanoiP[
#
〚i〛
,p-i+1]],{i,p-2}]+1}]&,t];First[Sort[v,Last[#2]>Last[#1]&]]]
In[]:=
superHanoi[{{d_},{a_,___,b_}}]:={d,a,b}​​superHanoi[{tower_,pegs_}]:=Module{a,pat,lp=Length[pegs],n,ans={},i,p,spread,back},​​a=Drop[hanoiP[Length[tower],lp],-1];​​pat=TableTaketower,1+
n-1
∑
i=1
a
〚i〛
,
n
∑
i=1
a
〚i〛
,{n,lp-2};​​spread=Table[p=Drop[pegs,{2,n}];i=Last[p];
p
〚-1〛
=
p
〚2〛
;
p
〚2〛
=i;{
pat
〚n〛
,p},{n,lp-2}];​​spread=Cases[spread,{{__},_}];​​back=Map[{First[#],Join[{Last[Last[#]]},Complement[Last[#],{Last[Last[#]],Last[pegs]}],{Last[pegs]}]}&,Reverse[spread]];​​Map[AppendTo[ans,superHanoi[#]]&,spread];​​AppendTo[ans,{Last[tower],First[pegs],Last[pegs]}];​​Map[AppendTo[ans,superHanoi[#]]&,back];​​Partition[Flatten[ans],3]

Utility- Animation code

Just put everything about the animation code altogether:

CITE THIS NOTEBOOK

Animate the optimal set of moves for 4-peg tower of Hanoi​
by Shenghui Yang​
Wolfram Community, STAFF PICKS, October 31, 2024
​https://community.wolfram.com/groups/-/m/t/3312046