Also, Silvia, nice work finding the WaveFunctionCollapse github, and the subsequent blog post. This is very well related to some of the other algorithms for growing Penrose tilings and to the Multicomputational Paradigm. Is anything like this available in Wolfram Function Repository?
One interesting property of the Penrose tilings (if I remember correctly) is: For any patch of a tiling also a toplogical disk, orientation and placement of all interior tiles is determined by the ring of tiles around the edge of the disk. The whole solution on the region is just a function of the matching rules, shape and chosen boundary condition.
In the case of non-tightly constrained tilesets, I would expect many different equivalent solutions with sufficiently large boundary. Sure the pictures turn out looking pretty, but it could also be an interesting combinatorial question how many patterns of square size
$n \times n$ can be formed from some chosen initial data?
The other way to go is try to make the algorithm more complicated, whether by machine learning / genetic algorithms, or by comparing with data. Maybe quantum chaos?