Given a set of rectangles $S := \{R_1, R_2, ..., R_n\}$, I want to place them on a canvas such that the bounding box contains as little unused space as possible. Rectangles must not overlap. Rotation is not allowed.
The canvas can be thought of as infinite in $x$ and $y$.
A rectangle $R_i$ is given by the coordinates of its left, top, right and bottom edges, $(l_i, t_i, r_i, b_i)$. The area a rectangle occupies is defined by all $(x, y)$ where $x >= l_i$ and $x < r_i$ and $y >= t_i$ and $y < b_i$.
To place rectangle $R_i$ on the canvas, choose a position $(x, y)$ and use this for the rectangle's top-left corner. That is, remove $R_i$ from the set, and add a new $R_i' := (x, y, x + (r_i - l_i), y + (b_i - t_i))$ to the set. The new rectangle $R_i'$ has the same width and height as the original $R_i$.
The bounding box is defined as a rectangle $B(S):=(l_{bb}, t_{bb}, r_{bb}, b_{bb})$ where $l_{bb}$ is the minimum of $l_i$ in all rectangles $R_i$ in $S$; $r_{bb}$ is the maximum of all $r_i$; $t_{bb}$ is the minimum of all $t_i$; and $b_{bb}$ is the maximum of all $b_i$. Thus, it depends on the positioning and size of each rectangle.
The unused space is the area of all $(x, y)$ inside $B$ for which there is no $R_i$ that overlaps that area. I want to miminise the unused space. Since the sizes of $R_i$ is fixed, minimising the unused space is the same as minimising the area of the bounding box $A(S) := (r_{bb} - l_{bb}) * (b_{bb} - t_{bb})$.
The algorithm should also prefer solutions whose bounding box resembles more a square than a rectangle. Thus, I want to minimise the score $Z(S) := k*A(S) + |(r_{bb} - l_{bb}) - (b_{bb} - t_{bb})|$. The score increases with the area of $B$ as well as with non-squareness of $B$ given by the difference of B's width and height. The factor $k$ is a constant that is large enough to make minimising the area much more important than the solution's squareness.
Although I'm interested in finding the optimal solution, this is not a requirement. I'm happy with an approximation, especially if it's simple.
(Application: I have a number of small texture images and I want to combine them to one large image that contains all of the given images, but does not use a lot more memory because of unused areas.)
One approach is to use mixed integer linear programming.
Your proposed objective function is nonlinear, but to get a linear objective function you could instead minimize the sum of width and height, with bounds on their difference to force the bounding box to be squarish.
See https://yetanothermathprogrammingconsultant.blogspot.com/2017/07/rectangles-no-overlap-constraints.html for details about how to impose constraints that prevent overlapping rectangles.