1. There is a SimplexNoise class in JOML. How can I use it to get transform values for the cubes?
You do not use noise to generate the transformation (i.e. "position") of a cube, but you use the noise function as a density function describing at every lattice point of your world whether or not that position contains a cube (and probably what kind of cube).
2. After getting the transform values, they must be stored somewhere. According to me it should be an array of model matrices. But Matrix4f returns a lot of errors and MatrixStackf does not provide Matrix array calculations.
You do not need matrices at all. They're overkill for what you are doing. Just assign each visible cube a 3-component vector denoting its position.
3. Now, how will I get my matrices to a buffer for shader use?
Like above, you do not need matrices. But if you wanted to, did you consider reading JOML's README.md
4. I need a simplified VAO/VBO initialization/binding/attrib code for the instanced arrays.
Google for "OpenGL instancing" and you'll find very good explanations and examples.