After some trial and error, here's what I think is the origin of the bug from my understanding. 
It's not a complete answer, but there's not enough space to write it as a comment.

The algorithm for the $p-$adic height follows, more or less, the algorith found in the paper by Mazur-Stein-Tate: (there are different conventions, e.g. on Sage the 1/p part is ignored etc, but the logic is the same)

[![enter image description here][1]][1]


The problem lies in the first step of the algorithm. In particular, we want to take $Q=mP$ such that:

(1) $Q$ reduces to $0\in E(\mathbb F_p)$ 

(2) $Q$ reduces to the non-singular component $E_{ns}(\mathbb F_\ell)$ over any bad prime $\ell$. 

To achieve the first, we have to multiply $P$ with $N_p= \# E(\mathbb F_p)$  (or if you like, with the order of $\tilde P$ in $\tilde E(\mathbb F_p)$ ). This part is fine, regardless the model.

For (2), if the elliptic curve is in global minimal model, it suffices to multiply $P$ with $c=\prod c_\ell$, so that $Q=cP$ will always reduce to a non-singular point (follows immediately from the definition...). When the Elliptic Curve is NOT in minimal model though, $Q=cP$ will not always reduce to a non-singular point. 

For reference, the reason we need (2) is so that  $$|d(nQ)|=|d(Q)^{n^2}f_n(Q) | \qquad (3) $$ holds for any n. (Here $f_n$ is the division polynomial and $d(Q)$ is the denominator as explained on step 2 of the above image). 

I am sure there are explicit ways of calculating the correct value for $c$ so that (2) (or equivalently (3) ) holds (I would be interested in knowing them).  A very  naive (but quick) way that I am currently using to resolve it is having a function that takes P=E.gen(0) and calculates the smallest value of c such that 
$$ |d(3 cP)| = |d(cP)^{3^2}f_3(cP) |$$

(I believe this is guaranteed to work as long as you take P as a **generator**. It is very inconvenient for complicated curves where finding generators is difficult/impossible, but for small curves the time is negligible. Hence a "correct" and explicit way of finding c would be helpful... I am not sure, mathematically, why or if that would imply (2) as I am only aware of the fact that $(2)\Rightarrow (3)$, not that a special case of (3) implies (2), but it seems that it does...) 

Defining $m$ as the least common multiple of $c$ and $N_p$ and taking $Q=mP$ should resolve the issue. I have tried it with a wide range of rank 1 elliptic curves and always get the expected value for the $p-$adic regulator that agrees both with the one from the minimal model and that also gives the correct prediction for lambda invariants.

This is not a complete answer as the described algorithm to find $c$ is a bit silly (and unproven) but I hope it provides some insight. Any corrections or improvements are welcome. 


  [1]: https://i.sstatic.net/o1c76.png