BUG: almost feasible for the restoration NLP #650
Answered
by
svigerske
TheRealelli
asked this question in
Q&A
-
|
Hello Everyone,
As far as I understand the message "Cannot call restoration phase at point that is almost feasible for the restoration NLP (violation 0,000000e+00)." is the import part. I have looked/googled for similar problems, but still couldn't figure it out. import java.util.Arrays;
import org.coinor.Ipopt;
public class Dice extends Ipopt {
// number of Periods
final static int T = 10;
final static int TLAST = T - 1;
final static int TFIRST = 0;
// Scalars adn Constants
final static double R = 0.03; // Rate of social time pref per year
final static double GL0 = 0.223;// Growth rate of population per decade
final static double DLAB = 0.195;// Decline rate of pop growth per dec
final static double DELTAM = 0.0833;// Removal rate carbon per decade
final static double GA0 = 0.15;// Initial growth rate for technology per dec
final static double DELA = 0.11;// Decline rate of technology per dec
final static double SIG0 = 0.519;// CO2-equiv-GNP ratio
final static double GSIGMA = -0.1168;// Growth of sigma per decade
final static double DK = 0.1;// Depreciation rate on capital per year
final static double GAMA = 0.25;// Capital elasticity in output
final static double M0 = 677;// CO2-equiv concentrations 1965 bill t C
final static double TL0 = 0.1;// Lower stratum temperature (C) 1965
final static double T0 = 0.2;// Atmospheric temperature (C) 1965
final static double ATRET = 0.64;// Marginal atmospheric retention rate
final static double Q0 = 8.519;// 1965 world gross output trillions 89 US dol
final static double LL0 = 3369;// 1965 world population millions
final static double K0 = 16.03;// 1965 value capital billions 1989 US dollars
final static double C1 = 0.226;// Coefficient for upper level
final static double LAM = 1.41;// Climate feedback factor
final static double C3 = 0.440;// Coefficient trans upper to lower stratum
final static double C4 = 0.02;// Coeff of transfer for lower level
final static double A0 = 0.00963;// Initial level of total factor productivity
final static double A1 = 0.0133;// Damage coeff for co2 doubling (frac GWP)
final static double B1 = 0.0686;// Intercept control cost function
final static double B2 = 2.887;// Exponent of control cost function
final static double PHIK = 140;// transversality coeff capital
final static double PHIM = -9;// Transversality coeff carbon ($ per ton)
final static double PHITE = -7000;// Transversalit coeff temper (bill $ per deg C
// PARAMETERS
// init parameters
static double[] l = new double[T]; // Level of population and labor
static double[] al = new double[T];// Level of Total factor productivity
static double[] sigma = new double[T];// Emissions-output ratio
static double[] rr = new double[T];// Discount factor
static double[] ga = new double[T];// Growth rate of T. F. P. from 0 to T
static double[] forcoth = new double[T];// Exogenous forcing other greenhouse gases
static double[] gl = new double[T];// Growth rate of labor 0 to T
static double[] gsig = new double[T];// Cumulative improvement of energy efficiency
static double[] dum = new double[T];// dummy variable 0 except 1 in last period;
// VARIABLES
String variableNames[] = { "MIU", "FORC", "TE", "TL", "M", "E", "C", "K", "CPC", "PCY", "I", "S", "RI", "TRANS",
"Y" };
// Setting Constants for easier access and readability
// MIU(T) Emission control rate GHGs 0
final static int MIU = 0;
// FORC(T) Radiative forcing, W per m2 1
final static int FORC = 1;
// TE(T) Temperature, atmosphere C 2
final static int TE = 2;
// TL(T) Temperature, lower ocean C 3
final static int TL = 3;
// M(T) CO2-equiv concentration bill t 4
final static int M = 4;
// E(T) CO2-equiv emissions bill t 5
final static int E = 5;
// C(T) Consumption trill US dollars 6
final static int C = 6;
// K(T) Capital stock trill US dollars 7
final static int K = 7;
// CPC(T) Per capita consumption thousands US dol 8
final static int CPC = 8;
// PCY(t) Per capita income thousands US dol 9
final static int PCY = 9;
// I(T) Investment trill US dollars 10
final static int I = 10;
// S(T) Savings rate as fraction of GWP 11
final static int S = 11;
// RI(T) Interest rate per annum 12
final static int RI = 12;
// TRANS(T) transversality variable last period 13
final static int TRANS = 13;
// Y(T) OUTPUT 14
final static int Y = 14;
// Setting Constants for easier access and readability
// EQUATIONS
// YY(T) Output 0
final static int YY = 0;
// CC(T) Cconsumption 1
final static int CC = 1;
// KK(T) Capital balance 2
// KK0(T) Initial condition of K 2 // not a real condition used for init KK(T)
final static int KK = 2, KK0 = KK;
// KC(T) Terminal condition of K 3
final static int KC = 3;
// CPCE(t) Per capita consumption 4
final static int CPCE = 4;
// PCYE(T) Per capita income equation 5
final static int PCYE = 5;
// EE(T) Emissions process 6
final static int EE = 6;
// SEQ(T) Savings rate equation 7
final static int SEQ = 7;
// RIEQ(T) Interest rate equation 8
final static int RIEQ = 8;
// FORCE(T) Radiative forcing equation 9
final static int FORCE = 9;
// MM(T) CO2 distribution equation 10
// MM0(T) Initial condition for M 10 // not a real condition used for init MM(T)
final static int MM = 10, MM0 = MM;
// TTE(T) Temperature-climate equation for atmosphere 11
// TTE0(T) Initial condition for atmospheric temp 11 // not a real condition used for init TTE(T)
final static int TTE = 11, TTE0 = TTE;
// TLE(T) Temperature-climate equation for lower oceans 12
// TLE0(T) Initial condition for lower ocean; 12 // not a real condition used for init TLE(T)
final static int TLE = 12, TLE0 = TLE;
// TRANSE(t) Transversality condition 13
final static int TRANSE = 13;
// Problem sizes
int variable_count = Y + 1; //Y last variable T
int n = variable_count * T;
int equation_count = TRANSE + 1;// Anz. individueller Bedingungen
int m = equation_count * T;
int n_j_one = 25;// short for nele_jac_for_period_one //some constraints are
// initialised with less variables
int n_j_more = 30; // nele_jac_for_more_then_one_period
int n_j_last = 6;// nele_jac_for_last_period
int nele_jac;
int nele_one;
int nele_hess;// from example
int count_bounds = 0;
int dcount_start = 0;
boolean printiterate;
public Dice() {
nele_jac = n_j_one + n_j_more * (T - 1) + n_j_last;
nele_hess = nele_one * T;
// calculat parameters in dependence of T
System.out.println();
for (int t = 0; t < T; t++) {
// GL(T) =(GL0/DLAB)*(1-exp(-DLAB*(ord(t)-1)));
gl[t] = (GL0 / DLAB) * (1 - Math.pow(Math.E, (-DLAB * t)));// because we start at 0 we add 1 to t
// L(T)=LL0*exp(GL(t));
l[t] = LL0 * Math.pow(Math.E, (gl[t]));
ga[t] = (GA0 / DELA) * (1 - Math.pow(Math.E, (-DELA * t)));
// AL(T) = a0*exp(GA(t));
al[t] = A0 * Math.pow(Math.E, (ga[t]));
// GSIG(T) = (GSIGMA/DELA)*(1-exp(-DELA*(ord(t)-1)));
gsig[t] = (GSIGMA / DELA) * (1 - Math.pow(Math.E, (-DELA * t)));
// SIGMA(T)=SIG0*exp(GSIG(t));
sigma[t] = SIG0 * Math.pow(Math.E, (gsig[t]));
if (t == TLAST) {
// DUM(T)=1$(ord(T) eq card(T));
dum[t] = 1;
} else {
dum[t] = 0;
}
rr[t] = Math.pow((1 + R), (10 * (1 - (t + 1))));// RR(T) = (1+R)**(10*(1-ord(t)));//because we start at 0
if (t < 14) {
// FORCOTH(T)$(ord(t) lt 15) = .2604+.125*ord(T)-.0034*ord(t)**2;
forcoth[t] = 0.2604 + 0.125 * (t + 1) - 0.0034 * Math.pow((t + 1), 2);
} else {
forcoth[t] = 1.42;// FORCOTH(T) = 1.42;
}
}
/* Index style for the irow/jcol elements */
int index_style = Ipopt.C_STYLE;
/* Whether to print iterate in intermediate_callback */
printiterate = false;
/* create the IpoptProblem */
create(this.n, this.m, this.nele_jac, this.nele_hess, index_style);
}
@Override
protected boolean get_bounds_info(int n, double[] x_L, double[] x_U, int m, double[] g_L, double[] g_U) {
// init bounds as no limit
for (int i = 0; i < n; ++i) {
x_L[i] = -2e19;// FORC,TL,CPC,PCY,S,RI,TRANS
x_U[i] = 2e19; // FORC,TL,M,E,C,K,CPC,PCY,I,S,RI,TRANS,Y
}
for (int i = 0; i < T; ++i) {
// * Upper and Lower Bounds:
// MIU.up(T) = 0.99;
x_U[i + MIU * T] = 0.99;
// MIU.lo(T) = 0.01;
x_L[i + MIU * T] = 0.01;
// K.lo(T) = 1;
x_L[i + K * T] = 1;
// TE.up(t) = 20;
x_U[i + TE * T] = 20;
// M.lo(T) = 600;
x_L[i + M * T] = 600;
// C.LO(T) = 2;
x_L[i + C * T] = 2;
// POSITIVE VARIABLES MIU, E, TE, M, Y, C, K, I;
// MIU,M,C,L can't be negative due to restraints above
// E
x_L[i + E * T] = 0;
// TE
x_L[i + TE * T] = 0 ;
// Y
x_L[i + Y * T] = 0 ;
// I
x_L[i + I * T] = 0;
}
// INIT
for (int i = 0; i < m; ++i) {
g_L[i] = 0;
g_U[i] = 0;
}
// all Equations equal to 0
// not true for KK(t) =L=, KC(TLAST) =L=, EE(T) =G= therefore:
for (int i = 0; i < T; ++i) {
// KK(t) //Except first period see after loop
g_L[i + KK * T] = -2e19;
// KC(TLAST)
g_L[i + KC * T] = -2e19;
// EE(T)
g_U[i + EE * T] = 2e19;
}
// Exept for the first KK
g_L[TFIRST + KK0 * T] = 0;
return true;
}
@Override
protected boolean eval_f(int n, double[] x, boolean new_x, double[] obj_value) {
assert n == this.n;
// UTILITY =E= SUM(T, 10 *RR(T)*L(T)*LOG(C(T)/L(T))/.55 +TRANS(T)*DUM(T));//GAMS is maximising
double utility = 0;
for (int t = 0; t < T; t++) {
utility += 10 * rr[t] * l[t] * Math.log(x[t + C * T] / l[t]) / 0.55 + (x[t + TRANS * T] * dum[t]);
}
obj_value[0] = utility;
return true;
}
@Override
// for constraint functions
protected boolean eval_g(int n, double[] x, boolean new_x, int m, double[] g) {
assert n == this.n;
assert m == this.m;
for (int i = 0; i < T; ++i) {
// YY(T) Output
// YY(T).. Y(T) =E= AL(T)*L(T)**(1-GAMA)*K(T)**GAMA*(1-B1*(MIU(T)**B2))/(1+(A1/9)*SQR(TE(T)));
g[i + YY * T] = x[i + Y * T] - (al[i] * Math.pow(l[i], (1 - GAMA)) * Math.pow(x[i + K * T], GAMA)
* (1 - B1 * (Math.pow(x[i + MIU * T], B2))) / (1 + (A1 / 9) * (Math.pow(x[i + TE * T], 0.5))));
// CC(T) Cconsumption
//CC(T).. C(T) =E= Y(T)-I(T);
g[i + CC * T] = x[i + C * T] - (x[i + Y * T] - x[i + I * T]);
// CPCE(t) Per capita consumption
//CPCE(T).. CPC(T) =e= C(T)*1000/L(T);
g[i + CPCE * T] = x[i + CPC * T] - (x[i + C * T] * 1000 / l[i]);
// PCYE(T) Per capita income equation
//PCYE(T).. PCY(T) =e= Y(T)*1000/L(T);
g[i + PCYE * T] = x[i + PCY * T] - (x[i + Y * T] * 1000 / l[i]);
// EE(T) Emissions process
//EE(T)..E(T)=G=10*SIGMA(T)*(1-MIU(T))*AL(T)*L(T)**(1-GAMA)*K(T)**GAMA;
g[i + EE * T] = x[i + E * T] - (10 * sigma[i] * (1 - x[i + MIU * T]) * al[i] * Math.pow(l[i], (1 - GAMA))
* Math.pow(x[i + K * T], GAMA));// k
// SEQ(T) Savings rate equation
//SEQ(T)..S(T) =e= I(T)/(.001+Y(T));
g[i + SEQ * T] = x[i + S * T] - (x[i + I * T] / (0.001 + x[i + Y * T]));
// RIEQ(T) Interest rate equation
//RIEQ(T).. RI(T) =E= GAMA*Y(T)/K(T)- (1-(1-DK)**10)/10 ;
g[i + RIEQ * T] = x[i + RI * T]
- ((GAMA * x[i + Y * T] / x[i + K * T]) - (1 - (Math.pow((1 - DK), 10))) / 10);
// FORCE(T) Radiative forcing equation
//FORCE(T).. FORC(T) =E= 4.1*(log(M(T)/590)/log(2))+FORCOTH(T);
g[i + FORCE * T] = x[i + FORC * T] - (4.1 * (Math.log(x[i + M * T] / 590) / Math.log(2)) + forcoth[i]);
if (i == TFIRST) {
// Initial Conditions
// KK0(T) Initial condition of K
// KK0(TFIRST).. K(TFIRST) =E= K0;
g[TFIRST + KK0 * T] = x[TFIRST + K * T] - K0;// TODO
// MM0(T) Initial condition for M
// MM0(TFIRST).. M(TFIRST) =E= M0;
g[TFIRST + MM0 * T] = x[TFIRST + M * T] - M0;
// TTE0(T) Initial condition for atmospheric temp
// TTE0(TFIRST).. TE(TFIRST) =E= T0;
g[TFIRST + TTE0 * T] = x[TFIRST + TE * T] - T0;
// TLE0(T) Initial condition for lower ocean;
// TLE0(TFIRST).. TL(TFIRST) =E= TL0;
g[TFIRST + TLE0 * T] = x[TFIRST + TL * T] - TL0;
} else if (i < TLAST) {
// KK(T) Capital balance
// K(T+1) =L= (1-DK)**10 *K(T)+10*I(T);
g[i + KK * T] = x[i + 1 + K * T] - ((Math.pow((1 - DK), 10)) * x[i + K * T] + 10 * x[i + I * T]);
}
if (i < TLAST) {
// MM(T) CO2 distribution equation
// MM(T+1).. M(T+1) =E= 590+ATRET*E(T)+(1 - DELTAM)*(M(T)-590);
g[i + 1 + MM * T] = x[i + 1 + M * T]
- (590 + ATRET * x[i + E * T] + (1 - DELTAM) * (x[i + M * T] - 590));
// TTE(T) Temperature-climate equation for atmosphere
// TTE(T+1).. TE(T+1) =E= TE(t)+C1*(FORC(t)-LAM*TE(t)-C3*(TE(t)-TL(t)));
g[i + 1 + TTE * T] = x[i + 1 + TE * T] - (x[i + TE * T]
+ C1 * (x[i + FORC * T] - LAM * x[i + TE * T] - C3 * (x[i + TE * T] - x[i + TL * T])));
// TLE(T) Temperature-climate equation for lower oceans
// TLE(T+1).. TL(T+1) =E= TL(T)+C4*(TE(T)-TL(T));
g[i + 1 + TLE * T] = x[i + 1 + TL * T] - (x[i + TL * T] + C4 * (x[i + TE * T] - x[i + TL * T]));
}
}
// KC(T) Terminal condition of K
// KC(TLAST).. R*K(TLAST) =L= I(TLAST);
g[TLAST + KC * T] = R * x[TLAST + K * T] - x[TLAST + I * T];
// TRANSE(t) Transversality condition
// TRANSE(TLAST).. TRANS(TLAST)=E=RR(TLAST)*(PHIK*K(TLAST)+PHIM*M(TLAST)+PHITE*TE(TLAST));
g[TLAST + TRANSE * T] = x[TLAST + TRANS * T]
- (rr[TLAST] * (PHIK * x[TLAST + K * T] + PHIM * x[TLAST + M * T] + PHITE * x[TLAST + TE * T]));
return true;
}
@Override
protected boolean eval_grad_f(int arg0, double[] arg1, boolean arg2, double[] arg3) {
// Auto-generated method stub
return false;
}
@Override
protected boolean eval_h(int arg0, double[] arg1, boolean arg2, double arg3, int arg4, double[] arg5, boolean arg6,
int arg7, int[] arg8, int[] arg9, double[] arg10) {
// Auto-generated method stub
return false;
}
@Override
protected boolean eval_jac_g(int n, double[] x, boolean new_x, int m, int nele_jac, int[] iRow, int[] jCol,
double[] values) {
assert (n == this.n);
assert (m == this.m);
int offset = 0;
if (values == null) {
for (int t = 0; t < T; t++) {
// YY(T) Output
// Y contains T, K, MIU, TE
// 4 entry
// sum of entrys = 4
iRow[offset] = t + YY * T;
jCol[offset++] = t + Y * T;
//
iRow[offset] = t + YY * T;
jCol[offset++] = t + K * T;
//
iRow[offset] = t + YY * T;
jCol[offset++] = t + MIU * T;
//
iRow[offset] = t + YY * T;
jCol[offset++] = t + TE * T;
//
//
// CC(T) Cconsumption
// CC(T).. C(T) =E= Y(T)-I(T);
// g[i+CC*T] = x[i + C * T] - (x[i + Y * T] - x[i + I * T]);
// CC contains, C, Y, I
// 3 entry
// sum of entrys = 4+3 = 7
iRow[offset] = t + CC * T;
jCol[offset++] = t + C * T;
//
iRow[offset] = t + CC * T;
jCol[offset++] = t + Y * T;
//
iRow[offset] = t + CC * T;
jCol[offset++] = t + I;
// CPCE(t) Per capita consumption
// CPCE(T).. CPC(T) =e= C(T)*1000/L(T);
// g[i+CPCE*T] = x[i + CPC * T] - (x[i + C * T] * 1000 / l[i]);
// CPCE contains: CPC, C
// 2 entry
// sum of entrys = 4+3+2 = 9
iRow[offset] = t + CPCE * T;
jCol[offset++] = t + CPC * T;
//
iRow[offset] = t + CPCE * T;
jCol[offset++] = t + C * T;
// PCYE(T) Per capita income equation
// PCYE(T).. PCY(T) =e= Y(T)*1000/L(T);
// g[i+PCYE*T] = x[i + PCY * T] - (x[i + Y * T] * 1000 / l[i]);
// PCYE contains: PCY, Y
// 2 entry
// sum of entrys = 4+3+2+2 = 11
iRow[offset] = t + PCYE * T;
jCol[offset++] = t + PCY * T;
//
iRow[offset] = t + PCYE * T;
jCol[offset++] = t + Y * T;
// EE(T) Emissions process
// EE(T)..E(T)=G=10*SIGMA(T)*(1-MIU(T))*AL(T)*L(T)**(1-GAMA)*K(T)**GAMA;
// g[i+EE*T] = 10 * sigma[i]
// * (1 - x[i + MIU * T] * al[i] * Math.pow(l[i], 1 - GAMA) * Math.pow(x[i + K * T], GAMA));
// EE contains: MIU, K
// 2 entry
// sum of entrys = 4+3+2+2+2 = 13
iRow[offset] = t + EE * T;
jCol[offset++] = MIU;
//
iRow[offset] = t + EE * T;
jCol[offset++] = t + K * T;
// SEQ(T) Savings rate equation
// SEQ(T)..S(T) =e= I(T)/(.001+Y(T));
// g[i+SEQ*T] = x[i + S * T] - (x[i + I * T] / (0.001 + x[i + Y * T]));
// SEQ contains: S, I, Y
// 3 entry
// sum of entrys = 4+3+2+2+2+3 = 16
iRow[offset] = t + SEQ * T;
jCol[offset++] = t + S * T;
//
iRow[offset] = t + SEQ * T;
jCol[offset++] = t + I * T;
//
iRow[offset] = t + SEQ * T;
jCol[offset++] = t + Y * T;
//
// RIEQ(T) Interest rate equation
// RIEQ(T).. RI(T) =E= GAMA*Y(T)/K(T)- (1-(1-DK)**10)/10 ;
// g[i+RIEQ*T] = x[i + RI * T] - (GAMA * x[i + Y * T] / x[i + K * T] - (1 - Math.pow((1 - DK), 10)) / 10);
// RIEQ contains: RI, Y, K
// 3 entry
// sum of entrys = 4+3+2+2+2+3+3 = 19
iRow[offset] = t + RIEQ * T;
jCol[offset++] = t + RI * T;
//
iRow[offset] = t + RIEQ * T;
jCol[offset++] = t + Y * T;
//
iRow[offset] = t + RIEQ * T;
jCol[offset++] = t + K * T;
//
// FORCE(T) Radiative forcing equation
// FORCE(T).. FORC(T) =E= 4.1*(log(M(T)/590)/log(2))+FORCOTH(T);
// g[i+FORCE*T] = x[i + FORC * T] - 4.1 * (Math.log10(x[i + M * T] / 590) / Math.log10(2)) + forcoth[i];
// FORCE contains FORC, M
// 2 entry
// sum of entrys = 4+3+2+2+2+3+3 = 21
iRow[offset] = t + FORCE * T;
jCol[offset++] = t + FORC * T;
//
iRow[offset] = t + FORCE * T;
jCol[offset++] = t + M * T;
if (t == 0) {// Initial COndition of Equations after the else
// KK0(T) Initial condition of K
// KK0(TFIRST).. K(TFIRST) =E= K0;
// g[INIT + KK0 * T] = x[INIT + K * T] - K0;
// KK0 contains: K
// 1 entry
// sum of entrys = 4+3+2+2+2+3+3+1 = 22
iRow[offset] = KK0 * T;
jCol[offset++] = K * T;
// MM0(T) Initial condition for M
// MM0(TFIRST).. M(TFIRST) =E= M0;
// g[INIT + MM0 * T] = x[0 + M * T] - M0;
// MM0 contains: M
// sum of entrys = 4+3+2+2+2+3+3+1+1 = 23
iRow[offset] = MM0 * T;
jCol[offset++] = M * T;
// TTE0(T) Initial condition for atmospheric temp
// TTE0(TFIRST).. TE(TFIRST) =E= T0;
// g[INIT + TTE0 * T] = x[0 + C * T] - T0;
// TTE0 contains C
// sum of entrys = 4+3+2+2+2+3+3+1+1+1 = 24
iRow[offset] = TTE0 * T;
jCol[offset++] = C * T;
// TLE0(T) Initial condition for lower ocean;
// TLE0(TFIRST).. TL(TFIRST) =E= TL0;
// g[INIT + TLE0 * T] = x[0 + TL * T] - TL0;
// TLE0 contains TL
// sum of entrys = 4+3+2+2+2+3+3+1+1+1+1 = 25
iRow[offset] = TLE0 * T;
jCol[offset++] = TL * T;
} else {
// KK(T) Capital balance
// K(T+1) =L= (1-DK)**10 *K(T)+10*I(T);
// g[i + KK * T] = x[i + K * T + 1]
// - ((Math.pow((1 - DK), 10)) * x[i - 1 + K * T] + 10 * x[i - 1 + I * T]);
// KK contains: K, I
// sum of entrys = 4+3+2+2+2+3+3+2 = 23
iRow[offset] = t + KK * T;
jCol[offset++] = t + K * T;
//
iRow[offset] = t + KK * T;
jCol[offset++] = t + I * T;
// MM(T) CO2 distribution equation
// MM(T+1).. M(T+1) =E= 590+ATRET*E(T)+(1 - DELTAM)*(M(T)-590);
// g[i + MM * T] = x[i + M * T]
// - (590 + ATRET * x[i - 1 + E * T] + (1 - DELTAM) * (x[i - 1 + M * T] - 590));
// MM contains: M, E
// sum of entrys = 4+3+2+2+2+3+3+2+2 = 25
iRow[offset] = t + MM * T;
jCol[offset++] = t + M * T;
//
iRow[offset] = t + MM * T;
jCol[offset++] = t + E * T;
// TTE(T) Temperature-climate equation for atmosphere
// TTE(T+1).. TE(T+1) =E= TE(t)+C1*(FORC(t)-LAM*TE(t)-C3*(TE(t)-TL(t)));
// g[i + TTE * T] = x[i + TE * T] - (x[i - 1 + TE * T] + C1 * (x[i - 1 + FORC * T]
// - LAM * x[i - 1 + TE * T] - C3 * (x[i - 1 + TE * T] - x[i - 1 + TL * T])));
// TTE contains: TE, FORC, TL
// sum of entrys = 4+3+2+2+2+3+3+2+2+3 = 28
iRow[offset] = t + TTE * T;
jCol[offset++] = t + TE * T;
//
iRow[offset] = t + TTE * T;
jCol[offset++] = t + FORC * T;
//
iRow[offset] = t + TTE * T;
jCol[offset++] = t + TL * T;
// TLE(T) Temperature-climate equation for lower oceans
// TLE(T+1).. TL(T+1) =E= TL(T)+C4*(TE(T)-TL(T));
// g[i + TLE * T] = x[i + TL * T] - (x[i - 1 + TL * T] + C4 * (x[i - 1 + TE * T] - x[i - 1 + TL * T]));
// TLE contains: TL, TE
// sum of entrys = 4+3+2+2+2+3+3+2+2+3+2 = 30
iRow[offset] = t + TLE * T;
jCol[offset++] = t + TL * T;
//
iRow[offset] = t + TLE * T;
jCol[offset++] = t + TE * T;
}
}
// KC(T) Terminal condition of K
// KC(TLAST).. R*K(TLAST) =L= I(TLAST);
// g[TLAST + KC * T] = R * x[TLAST + K * T] - x[TLAST + I * T];
// KC contains: K, I
iRow[offset] = TLAST + KC * T;
jCol[offset++] = TLAST + K * T;
//
iRow[offset] = TLAST + KC * T;
jCol[offset++] = TLAST + I * T;
// TRANSE(t) Transversality condition
// TRANSE(TLAST).. TRANS(TLAST)=E=RR(TLAST)*(PHIK*K(TLAST)+PHIM*M(TLAST)+PHITE*TE(TLAST));
// g[TLAST + TRANSE * T] = x[T + TRANS * T]
// - (rr[TLAST]) * (PHIK * x[T + K * T] + PHIM * x[T + M * T] + PHITE * x[T + TE * T]);
// TRANSE contains: TRANS, K, M, TE
iRow[offset] = TLAST + TRANSE * T;
jCol[offset++] = TLAST + TRANS * T;
//
iRow[offset] = TLAST + TRANSE * T;
jCol[offset++] = TLAST + K * T;
//
iRow[offset] = TLAST + TRANSE * T;
jCol[offset++] = TLAST + M * T;
//
iRow[offset] = TLAST + TRANSE * T;
jCol[offset++] = TLAST + TE * T;
}
return true;
}
@Override
protected boolean get_starting_point(int n, boolean init_x, double[] x, boolean init_z, double[] z_L, double[] z_U,
int m, boolean init_lambda, double[] lambda) {
assert init_z == false;
assert init_lambda = false;
if (init_x) {
// no lower bounds: FORC,TL,CPC,PCY,S,RI,TRANS
// no upper bounds: FORC,TL,M,E,C,K,CPC,PCY,I,S,RI,TRANS,Y
for (int i = 0; i < T; i++) {
x[i + MIU * T] = 0.02;// lower: 0.01 and upper: 0.99 //can't start at 0 due to bounds, arbitrarily set
// in bound
x[i + FORC * T] = 0;// lower: n.A. and upper: n.A.
x[i + TE * T] = T0;// lower: 0 and upper: 20
x[i + TL * T] = TL0;// lower: n.A. and upper: n.A.
x[i + M * T] = M0;// lower: 600 and upper: n.A.
x[i + E * T] = 0.01;// lower: 0 and upper: n.A. // set to 0.01, i wasn't sure if bounds are < or <=
x[i + C * T] = 2.01;// lower: 2 and upper: n.A. //can't start at 0 due to bounds, arbitrarily set in
// bound
x[i + K * T] = K0;// lower: 1 and upper: n.A
x[i + CPC * T] = 0;// lower: n.A. and upper: n.A.
x[i + PCY * T] = 0;// lower: n.A. and upper: n.A.
x[i + I * T] = 0.01;// lower: 0 and upper: n.A. // set to 0.01, i wasn't sure if bounds are < or <=
x[i + S * T] = 0;// lower: n.A. and upper: n.A.
x[i + RI * T] = 0;// lower: n.A. and upper: n.A.
x[i + TRANS * T] = 0;// lower: n.A. and upper: n.A.
x[i + Y * T] = Q0;// lower: 0 and upper: n.A. //TODO Q0?
}
}
return true;
}
protected void print(double[] x, String str) {
System.out.println("___");
System.out.println(str + ":");
for (int i = 0; i < x.length; ++i) {
System.out.print(str + "(" + i + ") " + ": " + x[i] + "; ");
}
System.out.println();
}
protected void printVariables(double[] x) {
System.out.println("VARIABLES");
for (int i = 0; i < variable_count; i++) {
System.out.println(variableNames[i] + i);
System.out.print(" ---beginn: " + (i * T));
System.out.print(" ---end: " + ((i * T) + T - 1));
print(Arrays.copyOfRange(x, (i * T), (i * T) + T), variableNames[i]);
System.out.println();
}
System.out.println();
}
/** Main function for running this code. */
// [main]
public static void main(String[] args) {
// Create the problem
Dice dice = new Dice();
// enable printing of current iterate in intermediate_callback
// dice.printiterate = true;
dice.setStringOption("hessian_approximation", "limited-memory");
dice.setStringOption("gradient_approximation", "finite-difference-values");
dice.setStringOption("jacobian_approximation", "finite-difference-values");
// Solve the problem
int status = dice.OptimizeNLP();
// Print the solution
if (status == SOLVE_SUCCEEDED) {
System.out.println("\n\n*** The problem solved!");
} else {
System.out.println("\n\n*** The problem was not solved successfully!");
}
double obj = dice.getObjectiveValue();
System.out.println("\nObjective Value = " + obj + "\n");
double variableValues[] = dice.getVariableValues();
dice.printVariables(variableValues);
double constraints[] = dice.getConstraintValues();
dice.print(constraints, "Constraint Values:");
double MLB[] = dice.getLowerBoundMultipliers();
dice.print(MLB, "Dual Multipliers for Variable Lower Bounds:");
double MUB[] = dice.getUpperBoundMultipliers();
dice.print(MUB, "Dual Multipliers for Variable Upper Bounds:");
double lam[] = dice.getConstraintMultipliers();
dice.print(lam, "Dual Multipliers for Constraints:");
}
}For anybody interested this is supposed to be the Dice 92-94 version translated from GAMS to Java IPOPT |
Beta Was this translation helpful? Give feedback.
Answered by
svigerske
Mar 15, 2023
Replies: 1 comment
-
|
If I enable option check_derivatives_for_naninf and increase the print_level to 7, I see something like this in the log: So the first equality constraint couldn't be evaluated and the current iterate. Further, the approximation of first derivatives is experimental and it is highly recommended to avoid using this. This usually gives very bad performance (many iterations or not solving the problem at all). With these two issues resolved, there is a much higher chance that Ipopt will converge. |
Beta Was this translation helpful? Give feedback.
0 replies
Answer selected by
svigerske
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
If I enable option check_derivatives_for_naninf and increase the print_level to 7, I see something like this in the log:
So the first equality constraint couldn't be evaluated and the current iterate.
This results in these
Cutting back alpha due to evaluation errorerror and you may want to reconsider your implementation to avoid these.For example, ensure that all functions can be evaluated and differentiated within the variable bounds, maybe dis…