Исключение NLopt nullptr в C при использовании NLOPT SLSQP (алгоритм на основе градиента)

Ребята из Servus, я делаю проект по идентификации параметров с помощью Nlopt(SLSQP), я написал тестовый код, но после 3 итераций компилятор всегда выдает исключение о Nullptr в 'grad' в 'myfunc'(objectfunction):

'Вызвано исключение: нарушение прав доступа для записи. град был nullptr.'

, где я использовал конечную разницу для вычисления градиента, потому что конечная разница может вычислить градиент сложной модели в моем проекте.

я попытался изменить другой шаг-размер с 1e-8 на 1e-6 для конечной разницы, после этого код работает нормально без исключения, я не знаю причину, может кто-нибудь сказать мне?

double myfunc(unsigned n, const double *x,double *grad,void *data){ 

    double h =1e-8;

    grad[0] = (log(x[0] + h) - log(x[0])) / h; //hier compiler throws exception
    grad[1] = (log(x[1] + h) - log(x[1])) / h;

    printf("\ngrad[0] is %10f grad[1] is %10f\n", grad[0], grad[1]);
    printf("\nx[0] is %10f x[1] is %10f\n",x[0],x[1]);

    return log(x[0]) + log(x[1]);
}

double myconstraint(unsigned n, const double *x, double *grad, void*data) {

    double *a = (double *)data; 
    grad[0] = a[0];
    grad[1] = a[1];
    return x[0] * a[0] + x[1] * a[1] - 5;
}

double myinconstraint(unsigned n, const double *x, double *grad, void *data) {

    grad[0] = 1;
    grad[1] = -1;
    return x[0] - x[1];
}

void main(){
    //test-code

    double f_max = -10000;
    double tol = 1e-16;
    double p[2] = { 1,2 };
    double x[2] = { 1,1 };
    double lb[2] = { 0,0 };
    double ub[2] = { 10000,10000 }; 

    nlopt_opt opter = nlopt_create(NLOPT_LD_SLSQP, 2);      
    nlopt_set_max_objective(opter, myfunc, NULL);

    nlopt_set_lower_bounds(opter, lb);
    nlopt_set_upper_bounds(opter, ub);
    nlopt_add_equality_constraint(opter, myconstraint, p, tol);
    nlopt_add_inequality_constraint(opter, myinconstraint, NULL, tol); 
    nlopt_set_xtol_rel(opter, tol);
    nlopt_set_ftol_abs(opter, tol);

    nlopt_result result = nlopt_optimize(opter, x, &f_max);//?

    printf("Maximum utility=%f, x=(%f,%f)\n", f_max, x[0], x[1]);

    system("pause");
}

hier - это результат в командном окне с размером шага 1e-8

град[0] равен 1.000000 град[1] равен 1.000000

x[0] is 1.000000 x[1] is 1.000000

град[0] равен 0,600000 град[1] равен 0,600000

x[0] is 1.666667 x[1] is 1.666667

град[0] равен 0,600000 град[1] равен 0,600000

x[0] is 1.666667 x[1] is 1.666667

после этого выдает исключение компилятора


person 王晓海    schedule 29.08.2019    source источник
comment
Каким должен быть третий аргумент nlopt_set_max_objective?   -  person Christian Gibbons    schedule 30.08.2019
comment
кроме того, main должен иметь возвращаемый тип int, а не void.   -  person Christian Gibbons    schedule 30.08.2019
comment
этот аргумент используется для импорта функции объекта и типа алгоритма (slsqp), max указывает Nlopt вычислить оптимальное значение параметра, когда функция объекта максимальна. кроме того, nlopt_set_min_objective может решить проблему минимальной оптимизации   -  person 王晓海    schedule 30.08.2019
comment
вы правы, я исправил это, но компилятор выдал то же самое исключение...   -  person 王晓海    schedule 30.08.2019
comment
log(x[0]) и double lb[2] = {0,0}; проверить, выдает ли журнал (0) ошибку   -  person J CHEN    schedule 30.08.2019


Ответы (1)


Вы должны проверить, является ли grad NULL, и вернуть градиент только тогда, когда он не равен NULL. Из документации:

Кроме того, если параметр grad не равен NULL, то мы устанавливаем grad[0] и grad[1] равными частным производным нашей цели по x[0] и x[1]. Градиент необходим только для алгоритмов, основанных на градиенте; если вы используете алгоритм оптимизации без производных, grad всегда будет NULL, и вам никогда не придется вычислять какие-либо производные.

Таким образом, ваш код должен выглядеть так:

double myfunc(unsigned n, const double *x,double *grad,void *data)
{ 
    double h = 1e-8;

    if (grad) {
        grad[0] = (log(x[0] + h) - log(x[0])) / h;
        grad[1] = (log(x[1] + h) - log(x[1])) / h;
    }

    return log(x[0]) + log(x[1]);
}
person user545424    schedule 27.09.2019