I think I have some suggestions.
You can make a binary variable $e_{ij}$ per edge corresponding to if those two vertices are in a clique together or not. Being in a clique is a transitive property. $e_{ij} \wedge e_{jk} = e_{ik}$. This "and" constraint is expressible as a MIP. This can be achieved with the constraint $e_{ik} <= e_{ij}$, $e_{ik} <= e_{jk}$, $
e_{ij} + e_{jk} - 1 <= e_{ik}$ link. This is the polytope generated by the truth table vertices (0,0,0) (0,1,0) (1,0,0) (1,1,1). This makes only vertex clique covers feasible. For maximizing edge weights you're good to go.
I don't have an elegant way to minimize the total number of cliques. Just maximizing the sum of all edges is a reasonable heuristic if heuristics are acceptable.
Another possibility is to take the sum of all the edges coming out of a vertex. This is the size of the clique it belongs to and is bounded by the total number of edges connected to that vertex.
$\sum_j e_{ij} = c_i - 1$, where $c_i$ is the size of the clique vertex $i$ belongs to. The sum over all vertices $\sum_j \frac{1}{c_i} = C$ is the total number of cliques $C$.
I don't have a good way to perform this inversion, except by MIP brute force.
This can be done by generating as many binary variables $a_{in}$ per vertex as it has neighbors, setting $\sum_n n a_{ni}=c_i$, enforcing that only one $a_{ni}$ is nonzero $\sum_n a_{ni} = 1$, and then making the objective $\sum_{jn}\frac{1}{n}a_{nj}$. This is rather dissatisfying to me and I hope there is a better way.
