How to solve QUBO problems in Q#?

To answer your first question, QAOA is an application of a hybrid classical-quantum algorithm, so using a classical optimizer is a perfectly valid solution here. I would suggest using e.g. scipy.optimize any other of your favorite optimization tools.

Your second question is regarding measuring the energy of the quantum state. Indeed, Q# does not allow inspection of the quantum state, because as you noted this is not possible on physical hardware either.

However, there are some ways to retrieve this info.

In an experiment on real hardware, what you would do is simply measure all qubits over and over again, let's say N times, which gives you a list of N bitstrings. If N is large this will give you a reliable probability distribution of all the qubit states. In simulation this can be done by measuring all the qubits in the qubit register in a for-loop, and counting the occurrences of each possible iteration. This should be straightforward to implement but I will suggest a solution below that uses existing library functionality which you might find helpful.

Instead of measuring the results of each qubit, we're going to implement a handy operation that lets us measure the probability of each valid state and the probability of obtaining an invalid state.

First of all, as explained in the blog post, there are the probabilities to get one of 3 valid states, and there's the probability of getting an invalid state:

// Allowed states
let state1 = [One, One, One, One, Zero, Zero];
let state2 = [Zero, One, Zero, One, One, One];
let state3 = [One, Zero, One, Zero, One, One];
let states = [state1, state2, state3];
mutable result = new Double[Length(states) + 1];

We'll keep track of the probability of getting an invalid state as the results come in by subtracting these values from 1.0.

mutable otherProb = 1.0; // Probability of any invalid state

So what we'll do is loop over all valid states, and then measure the probability of each, and return them and the probability of getting any invalid state in an array result.

for ((index, state) in Enumerated(states)) {
    let prob = MeasureProbabilityForState(state, numSegments, weights, couplings, timeX, timeZ, numMeasurements);
    set result w/= index <- prob;
    set otherProb -= prob;
}

set result w/= Length(states) <- otherProb;

The energy, or cost that we want to minimize, can then be calculated by multiplying the probabilities by the cost of each state.

(Obviously for this example this is a bit silly, as we can simply calculate the cost of each of the above allowed states and then pick the lowest value. However, in other quantum applications the space of valid states could for example be much larger, and in that case we would only need to evaluate the cost function for the states that have a nonzero probability.)

So, how to implement MeasureProbabilityForState?

We can use the operation EstimateFrequency for this. This operation measures the probability of getting a Zero result for given state prep and measurement operations. See docs here: https://docs.microsoft.com/en-us/qsharp/api/qsharp/microsoft.quantum.characterization.estimatefrequency.

The implementation can look something like this:

    operation MeasureProbabilityForState(
        state : Result[],
        numSegments: Int, 
        weights : Double[], 
        couplings : Double[], 
        timeX : Double[], 
        timeZ : Double[],
        numMeasurements: Int
    ) : Double {
        return EstimateFrequency(
            ApplyQAOA(_, numSegments, weights, couplings, timeX, timeZ), 
            MeasureRegisterIsInState(_, state),
            numSegments, 
            numMeasurements
        );
    }

where we have to check if the register is in the desired state and return Zero if it is:

    operation MeasureRegisterIsInState(register : Qubit[], state : Result[]) : Result {
        let result = MultiM(register);
        if (All<(Result, Result)>(EqualR, Zipped(result, state))) {
            return Zero;
        } else {
            return One;
        }
    }

and

    operation ApplyQAOA(
            x: Qubit[],
            numSegments: Int, 
            weights : Double[], 
            couplings : Double[], 
            timeX : Double[], 
            timeZ : Double[]
        ) : Unit {
        ApplyToEach(H, x); // prepare the uniform distribution
        for ((tz, tx) in Zipped(timeZ, timeX))
        {
            ApplyInstanceHamiltonian(numSegments, tz, weights, couplings, x); // do Exp(-i H_C tz)
            ApplyDriverHamiltonian(tx, x); // do Exp(-i H_0 tx)
        }
    }

This is just one way of using EstimateFrequency to solve this problem; I hope this example helps you figure out the best way that works for your quantum application.